Radial Domany-Kinzel Models with Mutation and Selection 
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We study the effect of spatial structure, genetic drift, mutation, and selective pressure on the 
evolutionary dynamics in a simplified model of asexual organisms colonizing a new territory. Under 
an appropriate coarse-graining, the evolutionary dynamics is related to the directed percolation 
processes that arise in voter models, the Domany-Kinzel (DK) model, contact process, etc. We 
explore the differences between linear (flat front) expansions and the much less familiar radial 
(curved front) range expansions. For the radial expansion, we develop a generalized, off-lattice DK 
model that minimizes otherwise persistent lattice artifacts. With both simulations and analytical 
techniques, we study the survival probability of advantageous mutants, the spatial correlations 
between domains of neutral strains, and the dynamics of populations with deleterious mutations. 
"Inflation" at the frontier leads to striking differences between radial and linear expansions. For a 
colony with initial radius Ro expanding at velocity v, significant genetic demixing, caused by local 
genetic drift, only occurs up to a finite time t* = Ro/v, after which portions of the colony become 
causally disconnected due to the inflating perimeter of the expanding front. As a result, the effect 
of a selective advantage is amplified relative to genetic drift, increasing the survival probability 
of advantageous mutants. Inflation also modifies the underlying directed percolation transition, 
introducing novel scaling functions and modifications similar to a finite size effect. Finally, we 
consider radial range expansions with deflating perimeters, as might arise from colonization initiated 
along the shores of an island. 
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I. INTRODUCTION 

To grow and prosper, populations must often migrate 
into new territory. These ubiquitous range expansions 
occur in bacterial growth on a Petri dish, tumor growth, 
viral infections, human migration, species movement in- 
duced by climate change, and in many other biological 
systems pHl]. Such expansions influence the evolution- 
ary dynamics of the population with, for example, en- 
hanced genetic drift due to low population densities at 
the frontier. To understand the universal features of 
these diverse range expansions, it is of interest to con- 
struct simple models with the essential elements of the 
population's evolutionary and spatial dynamics. 

Many features of a range expansion can be captured 
by adapting a simple "stepping stone" model d, M. In 
this model, individuals in a collection of island subpop- 
ulations, or demes, reproduce, die, and mutate stochas- 
tically. Two commonly used stochastic processes are the 
Moran and Wright-Fisher models (see PiJ-[l2| and refer- 
ences therein). The demes track spatial configurations, 
which simulate the spatial distribution of an evolving 
population. By allowing individuals to migrate between 
adjacent demes, one can simulate the effect of dispersal 
on the evolutionary dynamics in, say, one or two dimen- 
sions. Analysis of these models highlighted many impor- 
tant features distinguishing the evolutionary dynamics of 
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a well mixed versus spatially extended populations (see 
fTol | for a review). 

As discussed in Ref. [T3|, stepping stone models also 
provide a useful approximation to population genetics at 
the frontier of an expanding population. First, consider 
an asexual population expanding across a flat surface 
with an approximately linear front. Under a wide variety 
of conditions, the front will move with a constant speed 
v, set by the competition between growth and dispersal 
(see e.g., Ref. [u| and references therein). This growth 
is called a Fisher wave [l4| • Second, we assume that only 
individuals living near the population front are able to 
divide and settle into the new territory. In the approxi- 
mation of a flat population front, and in a reference frame 
moving with the frontier, we can simulate the front us- 
ing a one-dimensional stepping stone model. This model 
predicts that for small deme size, number fluctuations 
(genetic drift) create a strong spatial genetic demixing 
of the population into domains of closely related individ- 
uals. Initially, demixing occurs primarily through local 
fixation within a deme. At long times, however, coarsen- 
ing is dominated by the diffusive motion of the domain 
boundaries, which annihilate or coalesce on contact, cre- 
ating larger domains. Selective advantages will bias the 
motion of the boundaries, and mutations will introduce 
new domains into the population (Toj |. 

Stepping stone models (and the simplified voter-type 
models discussed below) of population genetics at fron- 
tiers make use of the notion of "dimensional reduction." 
A range expansion of organisms across a surface takes 
place in two spatial and one temporal dimensions. At 
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FIG. 1: A sketch of the phase diagram for mutation-selection 
balance in linear expansions with irreversible, deleterious mu- 
tations from green to red cells occurring at rate /x/, at a fitness 
cost s (see Appendix[S]for a biological motivation). The solid 
line between red and green regions represents the directed per- 
colation phase transition. We show as insets typical clusters 
generated from single green cells in the active and inactive 
phases. The dashed line at the top represents the case of no 
mutations where the system falls into the compact directed 
percolation (CDP) universality class. The dynamics at crit- 
icality (s = fij = at the black dot) in this CDP regime is 
shown. The dotted line shows the position of the mean field 
transition (i.e., for a well mixed population). An analogous 
plot for inflating radial range expansions is shown in Fig. 1111 



the cost of considerable computational complexity, a full 
2 + 1-dimcnsional model of the population genetics could 
be constructed, as was done recently within mean field 
theory to model range expansions of baker's yeast, Sac- 
charomyces cerevisiae [l5j]. However, the population den- 
sity and relative proportions of many organisms stop 
changing (and some even stop migrating) in the wake 
of a front colonizing virgin territory. The crucial dy- 
namics of genetic drift, competition, and spatial diffu- 
sion then takes place at a linear frontier, where one can 
approximate the 2 + 1-dimcnsional dynamics by a 1 + 1- 
dimensional stepping stone in a reference frame which 
moves with the population front. 

This mapping is only exact for perfectly flat population 
fronts. For neutral mutations, the flat front approxima- 
tion leads to boundaries of genetic domains behaving as 
simple random walks which wander from their initial po- 
sitions with a characteristic displacement Ax cx \ft that 
increases with the square root of time. Interesting com- 
plications arise when the front undulates - for neutral ge- 
netic variants, the boundaries between genetic domains 
couple to the undulations and wander more vigorously, 
which changes some of the detailed predictions of the the- 
ory [iol- ll6H18j . We neglect such complications here, but 
note that Kuhr et al. [19| have recently studied dynami- 
cal phase transitions in a two-species Eden model with an 
approximately flat but undulating front, and found dif- 
ferent critical exponents than those of the directed perco- 
lation universality class [20| considered here. It's worth 



noting that not only are flat one-dimensional habitats 
easier to analyze, a quasi-one-dimensional environment 
could describe the evolutionary dynamics along the bank 
of a river, along a seacoast, or along a constant altitude 
slope of a linear mountain range. To obtain inflation (or 
deflation) in this context, one could invoke, for example, 
a receding or advancing waterlinc around an island. Flat 
front models might also describe the gradual, climate- 
driven shift of some populations toward the Earth's poles 
II- 

An interesting special case of the stepping stone model 
is one in which a fit wild-type strain undergoes deleteri- 
ous mutations that confer a fixed selective disadvantage. 
If the deleterious mutations are irreversible, there will 
be a resultant selection-mutation balance, as shown in 
Fig. [T] The irreversible, deleterious mutations describe 
the genetic dynamics of a well adapted population evolv- 
ing in a fitness landscape with a single sharp peak (see 
Appendix [X] for more details) . The population sector 
dynamics exhibits a phase transition that falls into the 
directed percolation universality class [2(|. The order 
parameter is the fixation probability of the advantageous 
strain. In the "active" phase, the advantageous strain has 
a finite probability of surviving after an infinite amount 
of time (in a spatially infinite population) . In the "inac- 
tive" phase, the strain always dies out eventually. The 
transition into the inactive phase is sometimes called 
"mutational meltdown" in population genetics literature 

MM- 

In this paper we study curved population fronts in two 
dimensions that expand with some constant velocity v, 
motivated by populations that settle new territory from 
some localized, approximately circular, initial homeland 
with radius Rq. (A semicircular variant of this theory 
would describe organisms spreading into the interior after 
a coastal colonization event.) A particularly simple ex- 
ample is the growth of a circular bacterial inoculation on 
a Petri dish 22[ . We will again assume that the dynamics 
are confined to the population front. A dimensional re- 
duction to a 1 + 1-dimensional model will also be possible 
in this case, but the spatial dimension must now inflate 
to account for the growth of the curved front. To simplify 
the analysis, we assume the extreme limit of a stepping 
stone model, with a single cell per deme. As a result, the 
genetic drift, or number fluctuations, at the frontier will 
be very strong. We will argue that the long-time, large- 
size dynamics of this extreme limit of the stepping stone 
model is equivalent to the dynamics of a generalized voter 
model under an appropriate coarse-graining. 

In general, the curved front suppresses the effects of 
genetic drift by systematically increasing the perimeter 
of the frontier and inflating sector boundaries away from 
each other. The suppression becomes significant after a 
crossover time t* = Ro/v. A scaling argument for this 
crossover time follows from a comparison of the diffusive 
and inflationary length scales in the problem. At short 
times, sector boundaries move diffusively and sectors 
grow to a characteristic diffusive size id ~ -Ro^ ~ V Dt, 
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where (f> is the angular sector size, Rq is the homeland 
radius, and D is a diffusion constant. The diffusion con- 
stant scales as D ~ a 2 /r g where a is the cell diameter, 
and T g is the time between generations (see Sec. Mil for 
more details) . The inflationary growth due to the increas- 
ing population radius R(t) = Ro+vt has a characteristic 
size £i ~ 4>vt. So, there will be a crossover time t* at 
which the two length scales li and id are comparable. 
This sets id ~ Rq4> ~ \ /f Dt* ~ 0wt* ~ and we find the 
desired result i* ~ Rq/v. After time i*, sector bound- 
aries with a significant angular separation no longer dif- 
fusively interact, and the population separates into inde- 
pendent segments. Consequently, quantities such as the 
survival probability and the steady state fraction of an 
advantageous mutant become enhanced relative to linear 
expansions, because the advantageous mutant just has 
to survive until time t* after which inflation prevents its 
extinction via genetic drift. 

The inflationary dynamics also strongly influences the 
mutation-selection balance (directed percolation phase 
transition) by introducing a cutoff to the critical popula- 
tion dynamics after time t* . Just as in a finite size effect 
at a conventional phase transition, this cutoff is charac- 
terized by scaling functions. We find that the dynamics 
scales with t* in a similar way as for a finite popula- 
tion size L in a linear range expansion at criticality, but 
with a completely different scaling function. Moreover, 
inflation is qualitatively different from typical finite size 
effects due to transient dynamics or small system sizes 
[20I I23I ]: Inflation occurs at long times and is due to a 
gradually increasing system size. 

We will also consider the effects of inflation on the 
dynamics of single sectors of cells with a selective ad- 
vantage. A key dimensionless parameter that character- 
izes the strength of the selective advantage for voter-type 
models (in the presence of inflation) is K ~ sy / Ro/2vAt, 
where s < 1 is a selection coefficient, At is the time be- 
tween generations, and v is the circular front propagation 
speed. For k <C 1, selection is weak, genetic drift dom- 
inates the evolution, and the sector boundaries exhibit 
strong fluctuations. Conversely, when k 3> 1, the selec- 
tive advantage is strong, the sector boundaries move out 
deterministically at long times, and the fate of the sector 
is decided very early in the evolution. Examples of both 
cases are shown in Fig. [8] 

In general, if a subpopulation has a strong selective 
advantage, it will create "bulges" in an advancing pop- 
ulation front. The population front will also roughen 
due to the stochastic nature of the growth at the inter- 
face [lj| HH • These bulges and undulations will couple 
the population dynamics to the front geometry, changing 
the nature of the dynamics. As discussed above, these 
effects are not included in our study. However, we note 
that a smooth front approximation has been successfully 
employed in studies of expansions of bacteria, such as 
Pseudomonas aeruginosa In addition, this simplifica- 
tion will allow us to make connections between the radial 
expansions studied in this paper and the large body of 



work on stepping stone models for linear expansions with 
uniform fronts [lOj. 

The organization of the paper is as follows: In Sec. [TT] 
we develop a generalized Domany-Kinzel model [HJ to 
simulate both linear and radial range expansions. An un- 
fortunate artifact in lattice simulations of circular range 
expansions is genetic boundaries and colony shapes that 
lock into the discrete four- or six-fold symmetry of the 
underlying lattice. We describe novel simulation methods 
(which create an amorphous packing of organisms) that 
resolve this problem. In Sec. IIIII we interpret the simu- 
lation results by appealing to a generalized voter model 
[ID, description that includes both mutations and nat- 
ural selection. We also discuss how our generalized voter 
model is connected to the stepping stone model and the 
directed percolation universality class. The effect of in- 
flation on the directed percolation phase transition is also 
described in Sec. IIIII In addition, we compare linear and 
radial expansion results for quantities such as fixation 
probabilities and the heterozygosity correlation function. 
Sec. IIVI presents a brief discussion of an alternative "de- 
flating" inward expansion in which the population front 
gets smaller, as might happen in colonization from the 
perimeter of an island. We conclude with final observa- 
tions and a summary in Sec. [V] 



II. SIMULATIONS 

Many properties of the directed percolation phase tran- 
sition line shown in Fig. [T] have been studied using the 
Domany-Kinzel (DK) lattice model [HHIHI. The DK 
model is implemented on a hexagonal (i.e. equilateral 
triangle) lattice in which horizontal lines of sites are up- 
dated one at a time, starting from the top. The sites can 
be in one of two states (e.g., red or green) and the states 
of newly evolved sites are determined by the states of 
neighboring sites in the horizontal line above, as shown 
in Fig. HJa). Since all of the updating occurs at the bot- 
tommost actively evolving horizontal line, we have the di- 
mensional reduction discussed in the introduction. The 
downward direction on the lattice can be treated as a 
time axis. 

The DK model has a natural biological interpreta- 
tion. The sites represent cells, and the evolving lines 
of sites represent cell generations evolving as flat popula- 
tion fronts. Note that the DK model updates the entire 
front or generation of cells before moving the front. This 
scheme is called parallel updating [2(], |28[ and has non- 
overlapping generations. The two states of the cells in 
this model correspond to two genetic variants or "alle- 
les." During each update step, new "daughter" cells are 
evolved according to a probabilistic rule that depends 
on the states of neighboring "parent" cells in the previ- 
ous generation. Although these parents do not sexually 
reproduce, two neighboring asexual cells must compete 
for space for daughter cells during reproduction. Thus, 
not all of the parent cells will be able to propagate an 
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FIG. 2: (a) A sketch of the Domany-Kinzel (DK) model im- 
plemented on a hexagonal lattice so that each daughter cell 
has two potential parents. Periodic boundary conditions are 
employed, and each cell is evolved using input from the states 
of the two cells above it in the previous generation, (b) Il- 
lustration of the two update steps in our generalized DK 
model. The daughter shares the color of two parents of the 
same color. Otherwise, the daughter is green with probability 
PG = (2 — s) _1 ~ 1/2 + s/4 if its parents are different colors, 
where s < 1 is the selective advantage of the green strain. We 
set pn = 1 — pa — §E§ ~ 1/2 — s/4. A newly evolved cell 
is then allowed to mutate from green to red with a toward 
mutation probability (i/ < 1 or vice- versa with a backward 
mutation probability (ii<f 



radial range expansion 




FIG. 3: An illustration of a radial expansion simulation on a 
hexagonal lattice with an initial state given by a "homeland" 
of 7 cells. We add daughter cells one at a time, each time 
picking an unoccupied lattice site that is closest to the center 
of the homeland (the purple cell) . When two or more cells are 
at the same distance, we chose one of them at random. We 
implement the generalized Domany-Kinzel update rule that 
allows for more than two potential parent cells. Each daughter 
cell will be green with a probability pc that will depend on 
the number of green parents tig, given by Eq. |T]). Mutations 
can be implemented in a second update step as described in 
the text (also see Fig. [2]). 



associated daughter cell because of limited space at the 
frontier. As part of this reproductive process, successful 
daughters are allowed a small chance of switching colors 
via mutation, as discussed below. After a daughter cell is 
pushed forward into a frontier populated by a new gener- 
ation, cell division behind the advancing frontier ceases. 



A. A Generalized Domany-Kinzel Model 

In the original version of the DK model, there are two 
parameters p\ and p2 which set the update rules [24| ■ If a 
daughter cell has a green and red parent in the previous 
generation, the daughter is green with probability p\ . If 
both parents are green, the daughter is green with prob- 
ability p 2 - When the parents are both red, the daughter 
cell is always red so that no mutations from red to green 
cells are allowed. This parameterization is not well suited 
for population genetics, so a biologically motivated mod- 
ification that includes two-way mutations is necessary. 
To adapt and extend the DK model, daughter cells are 
now created in two steps shown in Fig[2j6). In the first 
step, a green or red daughter cell is created with a prob- 
ability pa or 1 — pa, respectively. We construct pc by 
drawing inspiration from the Moran model, a well studied 
update scheme for well mixed populations [IJ . The 
probability pg will be proportional to both the number 
of potential green parents uq , and the green cell growth 
rate, which is normalized to unity. The red cells suffer a 
selective disadvantage and grow at the smaller rate 1 — s. 



Then, to ensure that pq is properly normalized, we have 

riG n G 

PG = : c = : jz r, (1) 

nG + n R (l~ s) sn G + n T {l - s) 

where is the number of red parents and rif = uji+uq. 
After this first selection step, the surviving daughter cell 
mutates forward with probability \if <C 1 if it is green 
and backward with probability fib "C 1 if it is red. The 
update steps are shown for a simple linear expansion in 
Fig.Mfi), for which n T = 2. 

A lattice model for radial range expansions can now be 
constructed by taking advantage of the variability of tit 
in the update rule in Eq. (TTJ). For example, as illustrated 
in Fig. [3J we can use a hexagonal lattice and evolve the 
cells starting from an initial homeland centered around a 
reference cell. We always choose to evolve the daughter 
cell that is closest to the central reference cell in the initial 
seed. This choice ensures a uniform, circular population 
front. Certain daughters will now have three {tit = 3) 
instead of two (tit = 2) parents. Our revised DK update 
scheme easily adapts to this variation. As discussed be- 
low, this algorithm leads to circular colonies, and might 
arise biologically if previously deposited cells generate a 
chemical that stimulates cell division at the frontier. 



B. The Bennett Model 

Unfortunately, the radial DK model implemented on 
lattices with discrete rotational symmetries (e.g. hexag- 
onal, square, or triangular lattices) exhibits strong lat- 
tice artifacts. For example, domains of cells prefer to 
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grow along the crystallographic directions, which can in- 
troduce artificial oscillations in angular correlation func- 
tions, as discussed in Sec. IIIIBI To eliminate this prob- 
lem, we have implemented the DK model on an amor- 
phous lattice that does not have any special directions. 
An efficient way to generate such a lattice is provided by 
an extension of the Bennet algorithm (29j . 

In the original Bennet model [2|| , identical hard sphere 
particles (cells, in our case) were deposited sequentially in 
locations closest to a chosen reference point that becomes 
the center of the population. The possible locations for 
daughter cells are generated by enumerating all the pos- 
sible ways of placing a new cell so that it touches at least 
two cells without overlapping any other cell. The first 
few close-packed cells (a triangle of three, in the planar 
implementation discussed below) are placed manually in 
an initial "seed" . It turns out that, in two dimensions, 
the generated cluster of identical cells has a very strong 
hexa gon al ordering, regardless of the shape of the initial 
seed [30( | . Following Ref . [30| , we create an amorphous 
configuration by depositing cells with two different sizes 
with different probabilities. For a carefully chosen range 
of ratios p between the smaller and larger cell radius and 
the probability Q of depositing the smaller cell rather 
than the larger one, one finds packings with no preferred 
directions. 

One test of rotational symmetry for a generated cluster 
is to compute the structure factor 
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FIG. 4: (a) A radial DK model simulation on a hexagonal 
lattice. The structure factor of the lattice is shown in the 
right panel. The colors were assigned with neutral selective 
advantage (s = 0) and no mutations (pf — 0). The first 1500 
cells generated by the algorithm are shown. The shaded re- 
gion is the homeland of radius Ro ~ 10 cell diameters, where 
green and red cells are placed randomly with equal probabil- 
ity, (b) The same simulation performed on a Bennett model 
lattice with approximately equal numbers of small and large 
cells (Q — 0.5) and a radius ratio p = 0.769. The structure 
factor in the right panel shows that the Bennett call packing 
is isotropic. 



where N is the number of cells in the cluster and is 
the center of the i-th cell. A hexagonal cluster of cells 
with identical radii will have a striking structure factor 
with a clear six- fold symmetry, as shown in Fig. |4|a). 
By contrast, the amorphous cluster in Fig. 2^6) has a 
structure factor with no discrete symmetries. The latter 
cluster was generated using equal numbers of large and 
small cells (Q = 0.5) and size ratio p = 0.769. We use 
this choice of parameters in all radial range expansion 
simulations in the paper. 

Once an amorphous configuration of cells is generated, 
a homeland of cells is established by assigning genotypes 
of red or green to all cells within a radius Ro of the central 
reference point. Fig. 2J6) illustrates a well mixed initial 
homeland with Rq equal to about ten cell diameters. This 
initial homeland is then evolved using the rules discussed 
in Sec. Ill Al Similar to the evolution on the hexagonal 
lattice in Fig. [31 a new daughter is always the one that is 
closest to the central reference point used to generate the 
lattice. This can be regarded as an approximation for cell 
colonies where previous generations continuously secrete 
a chemical promoting cell division into the open space 
at the frontier. (We assume steric constraints prevent 
cell division in the interior of the colony.) The resulting 
concentration of stimulant will monotonically from the 
center, leading to circular colonies. Alternatively, a cir- 



cular population front can be imposed externally by, for 
example, a receding waterline around a circular island. 
The time t is measured in generations. Each genera- 
tion represents the time needed to evolve a new rim of 
cells at the circular population front. We assume that 
this rim has a finite thickness corresponding to the aver- 
age cell width, so that the population radius is given by 
R(t) = Rq + vt. The speed v is the rim thickness divided 
by the generation time. 

To facilitate the measurement of time in units of the 
radial colony size, we rescale all of our amorphous lat- 
tice coordinates such that we get the same cell number 
density as in a hexagonal lattice of identical cells with a 
diameter equal to one. We know that for the hexagonal 
lattice, a radial expansion will grow out by about a cell 
diameter per generation. Thus, with this rescaling, we 
can treat each generation in our amorphous lattice as all 
the cells in a rim with a thickness equal to one. This 
means that our population front speed is set to v = 1 
and the radius of the radial range expansion at time t is 
now R(t) = Ro+t. 

In the Bennett model lattice, there are a variety of 
ways to choose the competing parents in the local neigh- 
borhood of each daughter. In our simulations, we choose 
the parents to be the previously evolved cells that touch 



6 



the daughter cell. Since we update the cells one at a 
time, starting from the closest cell to the reference point, 
the adjacent parent cells must all be closer to the ref- 
erence point than the daughter. This usually gives us 
two or three potential parents. However, sometimes the 
daughter cell is adjacent to cells that are all further from 
the reference point and have not yet been updated. In 
these exceptional cases, we choose the competing parents 
to be the two closest nonadjacent cells that are closer to 
the reference point. 

Our amorphous cell packings will have slight inhomo- 
geneities. For example, there will be a variation in the 
coordination number z at each lattice site. The lattice 
spacing a will also vary across the lattice. Hence, we 
define an effective coordination number z r for our ra- 
dial expansions as well as an effective lattice constant a r . 
Since the amorphous lattice is not the densest possible 
packing of the cells, there will be gaps between the cells, 
and a r is expected to be slightly larger than the average 
cell diameter. In particular, we expect the amorphous 
lattice to be more loosely packed than a hexagonal lat- 
tice of identical cells with diameters equal to the average 
cell diameter in the amorphous lattice. 

The inhomogeneities, or "quenched disorder," in our 
amorphous cell packings require an average over spatial 
variations to get the effective parameters a r and z r . We 
have found it most straightforward to estimate a r and 
z r by fitting simulation results to known analytic solu- 
tions. For example, the heterozygosity function derived 
in Sec. MI Al is used to fit the combination a^/z r (see 
Fig. An alternative approach might be to use the 
coarse-graining procedure used to produce measured cor- 
relation functions from actual cell colonies in Ref. 
Some quantities, such as the fraction of green cells f(t) 
at time t, can be measured without knowing cither a r or 
z r . Note that the disorder is present in both the time-like 
(radial) and the space-like (circumferential) direction on 
the lattice. Fortunately, it can be shown [20j that spatio- 
tcmporally quenched disorder is an irrelevant perturba- 
tion in systems in the directed percolation universality 
class and it will not change the important features of the 
dynamics. 



III. VOTER MODELS WITH INFLATION 

To understand the radial simulation results analyti- 
cally, we shift to a generalized voter model description 
in which cells in the population are randomly updated 
one at a time, instead of generation by generation. The 
former and latter update schemes arc often referred to 
as random-sequential and parallel updates, respectively 
[28j . The original Domany-Kinzel model uses parallel 
updates as discussed previously. Our simulations (us- 
ing the DK model generalization discussed in Sec. HU 
represent an intermediate case since we evolve the sys- 
tem one cell at a time, but in a deterministic way to 
ensure that each generation is updated sequentially. The 



random-sequential scheme is used in canonical models 
such as the voter model or more general contact pro- 
cesses, which have features conducive to an analytical 
analysis [28j| . Parallel updates are easier to implement 
in simulations. Despite the differences between these 
two update schemes, we expect that they share universal 
features when we look at their coarse-grained dynamics. 
There are exceptional cases where this is not so, such 
as in the roughening transition in models of polynuclcar 
growth poj . However, we do not expect our simulations 
to be exceptional, because the directed percolation tran- 
sition can be modelled with either update scheme (20| . 
We start be considering flat frontier range expansions, 
and then introduce the inflation due to radial range ex- 
pansions in the coarse-grained description at the end. 

To start the analysis, we exploit dimensional reduction 
and consider the cells at the frontier of a (i+l-dimensional 
population as "voters" in a <i-dimensional space. It is 
convenient to describe the states of the cells with a set 
of spin variables {ci} ig z<i where <x; = ±1 and i is a d- 
dimensional vector of integers describing the location of 
the i th cell on a hypercubic lattice that approximates cell 
positions in our amorphous simulations. The spin vari- 
able encodes the cell state: er^ = I for a green cell and 
<7j = — 1 for a red one. No empty sites are allowed, which 
is equivalent to assuming a uniform front that advances 
in lock-step in (d + l)-dimcnsions with each generation. 
We now will use the notation {a} = {ci}igz d for conve- 
nience. The probability distribution of a particular spin 
configuration {cr}, P({a},t) obeys the master equation 

M 

d t P({a},t)=J2 Nab-K^Mi,*) 

-w {a}Mah P({o},t)], (3) 

where {cr}j is the same configuration of spins as {cr}, ex- 
cept that the j th spin is flipped (has an opposite state). 
Spin flipping in the biological context represents the re- 
placement of a cell at site j by a cell of the opposite 
type after cell division. The variable j in Eq. ^ runs 
over all spin locations and ui[ cr y^.[ tT y j are the probability 
rates of transitions between the configurations {<r} and 
{cr}j. The flipping rates co^y^^y . are given by the sum 
of three terms, 

, , , .drift I mut , scl ( a\ 

where Wr 1 "/' . , , is the rate due to random genetic 
drift, cjYy , i i is the contribution from mutations, and 
wf', , f , is the contribution from a selection bias. 

For neutral strains without mutations, only genetic 
drift contributes to the rates oj{ a y^ a y. . This genetic 
drift arises because a cell is more likely to change its state 
if it is surrounded by neighbors of the opposite type [3l[ , 
just as in the original voter model J2f|. Thus, uf^y^,-, . 



is the standard Glauber rate 1321 used in the voter model: 



drift 



2r„ 



i- u -i y , 

z ^ 
i n. n. j 



(5) 



where the summation is over the z = 2d nearest neigh- 
bors of voter (jj . The generation time r g corresponds to 
updating every cell once per generation. Mutations allow 
some cells to stochastically change their states indepen- 
dently of the states of their nearest neighbors. We model 
these processes by introducing the rates 



2r„ 



2r„ 



(6) 




FIG. 5: Three examples of surfaces generated by the set of 
points r(t,4>) given by Eq. (fl2|) with < t < 2 and — ty < 
4> < 7T. We set the initial population circumference and the 
crossover time to one: Rq = t* = 1. The = 1/2 surface 
is shaped like a bowl and has positive Gaussian curvature 
K > 0. Conversely, the = 2 surface is shaped like a trumpet 
and has K < 0. Confining a population to grow along the 
= 1 surface (with K = 0) is another way of implementing 
the regular radial range expansion with linear inflation. 



where \if and fib are the mutation probabilities for the 
cell during a generation time r g . These probabilities cor- 
respond to the mutation probabilities in the DK model 
simulations described in Sec. HU To make further contact 
with the DK model, green cells are replaced by red cells 
with a rate that is smaller than the reverse rate by a fac- 
tor of 1 — s, where s < 1 represents a selective advantage. 
The resultant contribution to the flipping rate is 



4-ZT, 



E 

i n . n . j 



(<Ti - 1) (1 + (Tj) . 



(7) 



where we sum over the nearest neighbors of spin ijj. 

We may now apply a coarse-graining procedure on the 
master equation (Eq. ©) with the rates in Eq. ^ to 
find a stochastic differential equation for the local, coarse- 
grained fraction /(x, t) of green cells: 



/(x,t) 



E 

ieBn(x) 



l + Oi(t) 



(8) 



where So is a set of spins in a neighborhood around fron- 
tier point x with volume ft. In Appendix[B]wc derive that 
for O = 1, /(x, t) obeys the equation 

cV(x, t) = D V 2 /(x, *) + 5 /(x, t) [1 - /(x, i)] 

- /2//(x, i) + a* 6 (1 - /(x, i)) + rjfx, t), (9) 

where Z) = a 2 /(zT g ) is the diffusion coefficient, s = s/t s , 
/2/ = fif /r g , jib = fJ-b/Tg, and 77 is a noise with correlations 
given by 



2/(x,t)[l-/(x,t)] 



xa^-iXx-x'). (10) 



This Langcvin equation can also be derived starting from 
a stepping stone model where the lattice sites are subpop- 
ulations undergoing a Moran process (see Appendix [3} 
and which exchange individuals between neighboring 
sites. By setting the number of individuals in each sub- 
population to one, we find the same Langcvin equation 



(Eq. ©), as shown in Ref. [T3|. Because of the non- 
linearity in the noise correlator Eq. (jTUJ) , the stochastic 
differential equation must be interpreted according to the 
ito calculus [13, HH . The time and space scaling used in 
the simulations sets r g = 1. 

The analysis so far is valid for both radial and linear 
range expansions as no reference has been made to the 
shape of the population front. For a linear front of length 
L in a two-dimensional population, the coordinate x is 
replaced by a single coordinate x € [— L/2,L/2]. For a 
circular front, we have to be more careful. In an "in- 
flating" front, the number N of active "voters" or cells 
around the perimeter increases in time. To cover a va- 
riety of interesting cases, we let the radius R(t) grow in 
time according to a general power law 



R{t) = Rq 



(11) 



where O is the growth exponent < < 00, t* is a 
characteristic time, and Rq is the initial homeland radius. 
In our simulations, we will be restricted to O = 1 where 
there is a constant front propagation velocity v such that 
t* = Ro/v. The choice of O = 1 is also appropriate for 
biological systems [3, [H, |33|, H3] expanding on a flat 
substrate. 

The 9 / 1 case might be used to model frontier ex- 
pansions on curved surfaces, such as a parabolic surface 
with positive Gaussian curvature, or a negatively curved 
surface like the horn of a trumpet as shown in Fig. [S] 
Each position r = (x, y, z) T on such a surface can be 
parameterized by the time t and an angle <j> £ [— it, tt). 
The angle <f> parameterizes a location on an expanding 
circumference. In Cartesian coordinates, 



(r 




Rq 



\ 



i + (f) 9 



COS 

sine/) 



(12) 



The population grows upward along the positive ^-axis, 
so that we can identify this direction with the time t. The 
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Gaussian curvature K(t, 4>) at each point (t, <f>) of such a 
surface can be calculated using standard methods [3f| : 

K(t, <f>) = * J ' . (13) 

i 2 +i?2(£) 20 82 

Note that K(t, <j>) is positive when O < 1 and negative 
when > 1. Surfaces at three different values of are 
shown in Fig. [5j 

The analysis will be general for now, and we will con- 
sider circular population fronts growing with arbitrary 
power 6 > of t, with radius R(t) given by Eq. (fTTT) . A 
location along the population front is given by the x and 
y components of r in Eq. (fT2|) : i?(i){cos <f>, sin </>}, where 
<j> G (— 7r, 7r) is the angle parameterizing a location on an 
expanding circumference. Setting x = R(t)(j> in Eqs. (|9]) 
and (|10j) . we obtain (upon setting = 0) the Langevin 
equation for a radial expansion 

W= *) = lJU 9 |/(^> *) + *) t 1 - /(^- *)] 



f{(j>,t)n f +Tl(<t>,t), 



(14) 



where the noise rj((f), t) obeys 



2a/(0,t)[l- /(<M)] 



5{t-t')6{ct>-(t>' 



(15) 

and we have set the generation time to one: r ff = 1. 
When we compare our simulation results to these voter 
models for radial range expansions, we will use an ef- 
fective lattice constant a r and coordination number z r , 
as discussed in Sec. HU The corresponding diffusion con- 
stant is then D = D r m a^/(z r T g ), with r g = 1. Similar 
equations to Eq. (|14j) and Eq. J 1 5 |) were derived via the 
stepping stone model in Ref. [l|, [lfj for the linearly in- 
flating case, where R(t) = Rq + vt. 



where ( )$> includes an angular average over <fi' and 
over realizations of the noise in the Langevin equation 
(Eq. ([Ti])) for /fai). The heterozygosity H(<p,t) repre- 
sents the probability of observing two cells of different 
types an angular distance </> apart at time t. Since our 
voter- type model has a single cell at each lattice point, 
-ff fa t) satisfies the boundary condition H (</> = 0, t) = 0. 

We set s = fj,f = in Eq. (fT4| and use the Ito calculus 
to get the equation obeyed by the angular heterozygos- 
ity H(cj),t) for the initial condition /(x, t = 0) = fo- 
We also assume that we can extend the range of <f> from 
the interval [— 7r,7r] to the entire real line (—00,00), an 
approximation will be valid as long as the average angu- 
lar sector size is small compared to 27r, the angular size 
of the population [3]. The heterozygosity H{4>,t) for 
an arbitrary power law inflation R(t) = Rq[1 + (t/t*) ] 
obeys a diffusion equation with a time-dependent diffu- 
sivity D(t) = 2D r /[R(t)} 2 . The time dependence can 
be removed by introducing a rescaled, "conformal" time 
variable 



t c (e,t) 



I 1 
e\i + (V**) e 

(6 - l) 2 Ji 



1 1 / 1 



(17) 



where 2-Pj fa P, 7i z ) is the Gaussian hypergeometric 
function [36l |. Changing variables to the conformal time 
yields a simple diffusion equation for H(<p,t c ), namely 



d tc H(^t c ) = 2D r dlH(^t c ). 



(18) 



Eq. (|18p can now be solved for the uniform initial condi- 
tion H{<j),t c = 0) = Hq = 2/o(l — fo) and the boundary 
condition H{(j> = 0,t c ) = 0. The solution, after changing 
back to the original time variable, is 



A. Neutral Evolution 

An important exactly soluble case is a neutral radial 
range expansion, considered previously for just linear in- 
flation (6 = 1 m Eq. ((H)) [J E3, El- We consider 
first the dynamics without mutations for simplicity. An 
illustration of a range expansion undergoing neutral dy- 
namics for a well mixed initial homeland with a fraction 
fo = 1/2 of green cells is given in the left panels of Fig.0] 
(a) and (b). We see in the figure that the domain walls 
between green and red sectors will initially annihilate in 
pairs due to genetic drift leading to domain coarsening. 
A useful measurement of the coarsening is the angular 
heterozygosity (for a spatially homogeneous initial con- 
dition) 



Hfai) = ^ + 0',t)[l-/^',t)] 



+ f{4>',t)[i-f{4> + 4>',t)] 



(16) 



Hfa t) = H erf 



Ro\4>\ 



y/8D r t c (&,t) 



(19) 



where t c (Q,t) is given by Eq. (fT7|) . 

We recover the known solution for a = 1 linear ex- 
pansion by replacing t c (0, t) with t in Eq. (fT9| [loj . The 
limiting value t c (Q,t — > 00) is particularly interesting as 
it tells us for which values of we can expect a "cutoff" 
in the dynamics, i.e., when the time variable t c (Q,t) ap- 
proaches a finite value and H fa t) approaches a limiting 
shape as t — > 00 in the radial expansion. In particular, 
we find 



t c (Q,t -> 00) 



DO < \ 

?r(e-i)r a 

2 sin(7r/0) 2 



(20) 



Thus, when the exponent is larger than the character- 
istic value 1/2 (which matches the diffusive wandering 
behavior of domain walls), the inflation will eventually 
dominate the diffusive dynamics. 
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FIG. 6: Nearly perfect collap se of the simulation data using 
the scaled variable f = Ro<j>/ y/9>D r t c (t) for the heterozygosity 
with Ho = 1/2, corresponding to random initial conditions, 
with equal probabilities for the two different colors. The ini- 
tial colony radius is given by Ro = 50, 100, 200, 300. Linear 
range expansions correspond to the limit Ro — > oo. All of the 
radial expansion simulations were fit with a single parameter 
D r . The scaling variable is £ = x/y/8Dt for linear expansions. 



The inflation-induced crossover also modifies the aver- 
age angular sector width A<j>(t). The width follows from 
Eq. (HU), upon noting that (see Eq. (C2) in Ref. [T(J) 



A0(t) 



f dH 



y/2irD r t c {t,e) 
RqHo 



(21) 



Eq. ([21]) will be valid only when A<p{t) <C 27T since we 
do not take into account sectors that take over the entire 
population. When 6 < 1/2, the angular width A<fi(t) 
diverges at long times so that only a single genetic vari- 
ant is present. Note the striking difference between this 
behavior and the > 1/2 case where the number of 
surviving sectors approaches a finite value given by 



A'siirv — 



2- 



lim. 



A<Xt) 



'2e 2 sin(7r/9) 
D r (@-l)t* ' 



(22) 



IK 



This result is in agreement with the results of Rcfs. 
03 for the case = 1 and t* = R /v. 

We now specialize to 6 = 1 and t* = Ro /v in order to 
check our results with simulations and known results. In 
this case, we have the simple relation 



t c (t) = t c (l,t) 



t 



t 



1 + t/t* 1 + vt/Ro 



(23) 



The transformation from t to t c (t) in Eq. (|23|) is a special 
case of the Mobius transformation and maps the interval 
t g (0,oo) to t c g (0,T) [13. Both the linear and radial 
expansion heterozygosities collapse to a single universal 
function H(£) = i? erf [^| 5 with £ = x/y/SDi for lin- 
ear expansions, and £ = Ro<p/ ^%D r t c (t) for radial ones. 



These theoretical predictions for the heterozygosity are 
described in more detail in Ref. [l(J and are compared to 
simulation data in Fig. [6l The comparison requires fit- 
ting the radial diffusion constant D r = a%/(z r T g ), which 
depends on the effective radial lattice constant a r and 
the coordination number z r for the Bennett model cellu- 
lar packing (r s = 1 in the simulations). We find excellent 
data collapse and agreement between simulation and the 
voter model theory. 



B. Heterozygosity in Radial Range Expansions 

It is also possible to derive analytic solutions for the 
heterozygosity in the biologically relevant limit of neutral 
mutations. Neutral mutations are common and can have 
an important effect on a population's genetic composi- 
tion [381 ] . We now compare strictly neutral mutational 
dynamics with genetic drift in radial versus linear range 
expansions. Man y p revious results for the linear case are 
reviewed in Ref. [lCf . 

Mutations spoil the clean connection between linear 
and radial range expansions via a conformal time coor- 
dinate (Eq. (|23[) h Instead, mutations (we now assume 
both [if and \x b are nonzero) dominate the dynamics at 
long times in both types of range expansion. The equa- 
tion of motion for H(<p, t) can be derived from Eq. (fT4|) 
via the Ito calculus, just as in the neutral case. An ad- 
ditional term fj,bf(x,t) is added to the left-hand side of 
Eq. (|14[) to take into account backward mutations from 
red to green cells at rate fib- The equation of motion for 
a radial range expansion for a homogeneous initial cell 
density /(x, t = 0) = f is 



d t H{(j>,t) 



2D 

——dlH(ct>, t) - 2(j* f + [i b )H{^ t) 



+ 2(fi f -fi b )f(t) + 2fi b , 



(24) 



where the third and fourth terms on the right-hand side 
of Eq. [24] come from the coupling between the equations 
for the first and second moments of /(x, t). These terms 
include a contribution from the average fraction of green 
cells /(<) = (/(x,t)) x , given by 

m = ^ + [/o(M/+M6)-He- t( ^ + ^ ) i (25) 
M/ + M6 

(see Ref. [T3|) which decays to a steady state value f(t — > 
00) = /it,/(/i/ + /ifc). To get the analogous linear equation, 
one can replace the angular coordinate 4> with a spatial 
coordinate x along the linear front and set R(t) — >■ 1. 

The solution to Eq. ([2~4"[) for an initially homogeneous 
random initial condition H(4>, t = 0) = Hq = 2/o(l — jo) 
is 

H((f>,t) = H^KRo + vt)cj)} + e- { ^ +t * b)t H tr {cj),t), (26) 
which at long times approaches the steady state 
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#oo[(-Ro + vt)(f>], where 



l_ e -WV^ 



(M/ + Mb) 2 
The transient term Ht v (4>,t) is given by 



(27) 



fftrfa,*) =H e- t ' erf 



yi + t/**_ 

M/ + Mb 



x / erf 
'o 



y/(l + t*/t){l-y) 



-t' /■! 



A*/ + Mb Jo 
2e*' fiffib 
(M/ + Mb) 2 



erf 



^/(l + t*/t)(l-y) 



e^dy 
e 2t ' y dy 

(28) 



where we have defined a dimensionless time i' = (/x/ 
/ib)i, and a dimensionless length scale 



(go + ^ 



(29) 



The steady-state shape -ffoo[(-Ro + «*)</>] is exactly the 
same for linear range expansions with mutations (see 
Ref. [l(|) if we identify (i? + u*)0 = x as the separa- 
tion between two cells in the linear expansion. In both 
expansions, mutations set the average domain size at 
long times. Also, a careful asymptotic analysis of the 
integrals in Eq. (|2"8")l reveals that the transient term in 
Eq. (f2"o| decays as e~ < - Pf+tJ " b ) t . A similar analysis shows 
that Eq. (|28|l reduces to the correct transient term for a 
linear range expansion in the limit t* — ¥ oo. More de- 
tails about the linear range expansion heterozygosity are 
discussed in Ref. [l(J . 

Eq. (|2"B"|) reveals two competing time scales: the mu- 
tational relaxation time (/iy + A*b) 1 and the inflation 
time t* — Ro/v. Mutations equilibrate on the first time 
scale, while the second determines when inflation sets in. 
The effect of inflation is observable for intermediate time 
scales such that t* <C t <C (/i/ + Mb) -1 - The average 
angular domain size A.<f>(t) = [c^i? (</>, t)|^ =0 ] _1 in this 
intermediate time regime is given by 



A0(i) 



HqRq 



1 - 2t' 



+ 



7T 
^0 



Mb (2/p - 1) 
Mb + M/ 




0(r/ 2 ) 



(30) 



where i' = (m/ + Mb)* ^ 1 an d *c(*) is the conformal 
time coordinate (Eq. ([22])). The factor in front of the 
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FIG. 7: A calculation of the heterozygosity H{cj>,i) at time 
t = 60 2> t* on different lattices. i?o = 15, v = 1, /if = 0.1, 
and s — 0.57 in the simulations. The length and time units are 
measured using the rescaling discussed in Sec. HU Oscillations 
in the heterozygosity on regular lattices reflect the 4— and 
6— fold rotational symmetries. These oscillations are absent 
for the amorphous Bennett model lattice with p — 0.769 and 
Q = 0.59. 



curly braces in Eq. ([50)1 is the neutral result of Eq. (f2"Tj) 
with = 1. 

The heterozygosity correlations are sensitive to dynam- 
ics of our model on different lattices. When we are in the 
active phase of the dynamics with selection and one-way 
mutations (/ij, = 0) on a regular lattice, the heterozygos- 
ity function exhibits oscillations as shown in Fig. [7J The 
number of oscillations reflects the rotational symmetry 
of the lattice, as shown for the hexagonal (6 oscillations) 
and square (4 oscillations) lattices in Fig. [7l The oscilla- 
tions arise because the green cell domains grow preferen- 
tially along the crystallographic directions of the lattice. 
These lattice artifacts obscure the evolutionary dynam- 
ics. As discussed in Sec. HH the problem is corrected by 
using an amorphous lattice which smooths out the oscil- 
lations. 



C. Single Seed Dynamics and Selective Advantage 

Finding exact solutions is difficult when we include se- 
lection. Unlike the cases considered in the two previous 
sections, it is not possible to apply the Ito calculus to 
Eq. (|l~4l to exactly calculate the angular heterozygos- 
ity H(<j>,t). The nonlinear term s/(x, t)[l — /(x, t)] in 
Eq. (|14[) leads to an infinite set of equations involving 
ever higher order moments of the field /(x, t). Although 
various approximations can be tried, accurate analytic 
results are limited [l(j. Hence, to study the effects of 
selection on radial range expansions we instead consider 
the fate of a "seed," i.e. a single initial green sector of 
size 4>o <C 1 with selective advantage s. We assume that 
mutations are rare, so that /jj = ^ « 0. During each 
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and 



FIG. 8: Two examples of single sector dynamics with Ro = 50 
(in units of the average cell diameter) evolved over about 250 
generations for two very different values of the key inflationary 
parameter k ~ s^/R Q /2vT g discussed in the text. In (a) a 
combination of selection bias and inflation dominates genetic 
drift after a short initial transient, with k ~ 2.4. The initial 
sector width in radians at the edge of the homeland is <j>o w 
0.1. In (b) Hi ~ 0.05 <C 1, so genetic drift can influence the 
result expected from selection bias and inflation alone for an 
extended period. We set <j>o ~ 0.2 in this case. We see that 
the fate of the sector in (a) is decided much faster than the 
fate of the sector in (&), which still might eventually die out. 



generation (corresponding to a time step r s ), the size <f> 
only changes due to the competition between red and 
green cells at the two boundaries of the sector. We make 
the further simplifying assumption, easily imposed in our 
simulations, of neglecting the bulge expected for a seed 
representing a favorable mutation at the frontier [l5|, [l8| . 
This approximation should be adequate for < s < 1, 
and it allows analytic progress for nearly circular colonies. 

Each generation in the radial expansion simulations 
consists of a rim of cells with an approximately single 
cell thickness (see Sec. |TTJ) . Hence, the green sector size 
changes each generation due to a single green cell com- 
peting with a neighboring red cell at each of two bound- 
aries. Specifically, each boundary can move by some av- 
erage amount ±A(f> w ±a r /R(t) per generation, where 
a r is an effective spacing between adjacent cells along 
the circular population front approximately equal to the 
effective lattice spacing a r (see Appendix [Cjl. The prob- 
ability distribution for the angular size <j> of a sector sat- 
isfies the Fokker-Planck-like equation (also known as a 
Kolmogorov forward equation - see Ref . [Io[ ) 



d T Ps 



w 



1 



dp s 
d(f> 



d 2 p s 
dtp 



(31) 



where t = t c (t)/t* = vt/(Ro + vt) = vt/R(t) is a dimen- 
sionless conformal time coordinate. The term in Eq. (|3"Tj) 
proportional to w represents a bias due to selection, and 
the term proportional to A represents genetic drift. As 
shown in Appendix [C] w is proportional to the selective 
advantages enjoyed by the green cells, 



(for s < 1), 



(32) 



A 



RnVT a 



(S « 1), 



(33) 



where 7 is a lattice-dependent factor of order unity. For 
our disordered Bennett model lattice, we find a r w 0.90 
(in units of the average cell diameter) and 7 ~ 1.2. We 
impose the absorbing boundary condition p s (4> = 0,r) = 
0, just as in the linear case [13]. Note from Eq. ([31"j) 
that inflation suppresses both drift and diffusion by fold- 
ing the entire time evolution of the range expansion into 
the conformal time r. Inflation suppresses the diffusion 
more than the drift term since they scale with R(t)~~ 2 
and R(t)~ l , respectively. 

The Langevin equation corresponding to Eq. (|31l) is 



dr 



(34) 



where r](t) is a Gaussian white noise (with (r](t)r)(t?)) = 
5(t — t')). Note that unlike the multiplicative noise in 
the Langevin equation for the coarse-grained cell density 
/(x, t) (Eq. ©), Eq. (|3"4"|) has a simple noise term. Hence, 
single sector dynamics is more easily analyzed using the 
"dual" description of sector boundary motion, instead 
of the full density /(x, i) dynamics in Eq. ([§]) (also see 
Appendix ICj) . When R ^> a r , we might hope to neglect 
genetic drift relative to the bias term. Upon setting A = 
0, Eq. §5Q is solved by 



4>{t) = 00 



w In (1 — t 
-R(t) 



< In 



Ro 



(35) 



similar to the logarithmic spiral sector boundaries found 
in Ref. [l5[ . To obtain the conformal time to fixation 77 , 
we set 4>( T f) = 27r to find 



Tf = 1 — exp 



2tt 



(36) 



for the purely deterministic model. When noise is in- 
cluded in our radial voter model, we find from our simula- 
tions that Tf w 1 for all s < 0.5 and cf>o < it. Thus, a sec- 
tor with small initial angular width 0o <C 1 can only take 
over a population in a reasonable time provided ai « s 
is large. When <fio is very small, the takeover predicted 
by the deterministic solution in Eq. ([35]) is preempted by 
number fluctuations, i.e., genetic drift. If a small initial 
sector survives the genetic drift, it can eventually take 
over the entire population in a deterministic fashion. 

To go beyond the deterministic approximation, we 
must account for the diffusive dynamics of the sector 
boundaries caused by genetic drift. To treat the case 
A > 0, we define a new angle variable 



ip{r) = 4>{t) + wln(l — t), 



(37) 
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FIG. 9: A plot of the crossover time r c that marks when 
genetic drift is no longer able to keep up with inflation. The 
time t c is given as a dimensionless time coordinate, related to 
the time t via r = vt/(Ro + vt) = vt/R(t). The dashed lines 
show the asymptotic approximations to r c given in Eq. flTj) . 



such that the sector probability distribution p s (ip, t) now 
obeys 



d T p s (ip,T) = A 



d 2 p s {jj 1 T) 
dip 2 



(38) 



Our boundary condition p s {4> = 0, r) = moves in time, 
p s {ip = wln(l — t),t) = 0. The diffusive motion of the 
variable ip, with uncertainty 6ij) = ±v2Ar, will be un- 
able to keep up with the absorbing boundary condition 
at ip = w ln(l — r) after a crossover time r c , such that 



w 



V2A 

with the key parameter 



ln(l - t c ) = -Kln(l - t c ), 




(39) 



(40) 



A plot of t c (k) is shown in Fig. [9l The behavior in the 
two limits k> 1 and n -C 1 is given by 



7"c(«) 



1 - e~ x l K k < 1 

K~ 2 K > 1 



(41) 



When k ^> 1 (strong selection), the transition from 
trajectories dominated by genetic drift to those domi- 
nated by deterministic natural selection is very fast. The 
inflationary effects on the survival probability are then 
negligible, and the relevant time scale is the diffusion 
time Ti ~ </>q/(2A). For times r 3> t\, the survival prob- 
ability of a sector for radial voter models is described by 
the long-time linear range expansion result [l8j 
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FIG. 10: Radial survival probabilities Sr(r) for a green sector 
(like that in Fig. [8]) of size cf>o ~ 0.036 initialized on the rim of 
a homeland with Ro = 300, evolved over 1750 generations, for 
a variety of selective advantages s. The solid lines represent 
the adiabatic approximation to S r (r) (Eq. (|46|) V The survival 
probabilities were averaged over 8 x 10 4 independent runs. 
We fit the black line with the lattice parameter a r , finding 
a r m 0.9. The other lines use 7 w 1.2 (see Appendix [C)) . 



If r <C t 1 , however, we must either use the full linear 
survival probability solution of Eq. (f3"Tj) (Eq. 3.2.15 of 
Ref. Hl|), or else exploit the small time solution 



S r {T < n) w 1 



4Ar 



■cxp - ' - - - — - (k>1). 

(43) 

In our simulations, it is difficult to access the regime 
k ^> 1 for small s since this requires large values of the 
initial radius Rq. In practice, we arc limited to the range 
< k < 10. Fig. E[a) shows a simulation for k ~ 2.4, 
with only minor diffusion motion in the sector boundaries 
visible at the onset of the range expansion. In this case, 
if the sector survives beyond the relatively small time r c , 
it is quite likely to survive indefinitely. 

When k <C 1, however, genetic drift dominates the 
deterministic sector size prediction in Eq. (|35p over the 
entire interval t € (0, r c ). As discussed above, a single 
sector is quite unlikely to take over the whole population 
for values of </>o not too close to 2-7T. We can then approxi- 
mate the diffusive motion of ijj = 4>+w ln(l — t) by confin- 
ing it to a semi-infinite interval between b(r) = w ln(l— r) 
and co, instead of imposing an upper limit of 27r. 

When r < t c , the boundary condition p s (ip — w ln(l — 
t),t) = moves slowly compared to the diffusion, al- 
lowing an adiabatic approximation (26j : The Fokker- 
Planck equation (Eq. (pTTjO for the sector size probability 
distribution p 3 (<f>,T) is treated as a drift-diffusion equa- 
tion with diffusion coefficient A and a time-dependent 
drift velocity \db(T)/dr\ = w/(l — r), evaluated instan- 
taneously at r. With the initial probability distribution 
p s (4>, 0) = 6((f> — 4>o), the solution to Eq. (pTTj) in this ap- 
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proximation follows from the method of images [26|: 
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(44) 



The change dS r (r) / dr in the radial survival probability 
must be given by the probability flux [Ad<^p s (0, r)] |^_ 
evaluated at the absorbing boundary. Integrating this 
probability flux for a radial expansion yields 



S r (r) = l-A / dr' [3^ s (0,t') 
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(46) 



where an asymptotic (small A and u>) approximation is 
used to evaluate the integral in Eq. (|45[) . Fig. [TU] shows 
that Eq. (|46[) is a remarkably good approximation to 
S r (r) even in the selection dominated regime r > r c . 

In the weak selection limit re — > 0, the movement of 
the boundary condition can be ignored entirely. Eq. (|46[) 
then reduces to the survival probability of a diffusing par- 
ticle on a semi-infinite interval with an absorbing bound- 
ary at the origin [26[ : 
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Fig. [TU] shows that the theory matches the simulation re- 
sults for radial expansions with O «a 0.036 and R = 300 
well, especially at small values r = vt/(Ro + vt). As 
t — > 1 (i.e. in the limit of long times), both inflation 
and selection prevent the extinction of the sector and 
S r (r) approaches a nonzero limit for all values of s, in- 
cluding s = 0. An asymptotic analysis of the integral in 
Eq. (|45[) for small ai a s shows that the limiting value 
S r (r 1) = Soo for small initial angular sector size </>o 
is approximately 
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For completeness, we also consider the regime = 
27r 2 /A <C 1, where genetic drift is strong enough to al- 
low an initially small green sector in Fig. [8] to diffusively 
take over the entire population. In this case, the finite 
range of the sector size <p becomes important. In the ab- 
sence of selection (re — > 0), the survival probability can be 
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FIG. 11: The color map of the average fraction of green 
cells at long times (t = 1800 generations) for a radial ex- 
pansion with a homeland radius of Ro = 200, averaged over 
400 runs, for positive values of pf and s. The average frac- 
tion p(t) approaches a steady state at long times so that 
Poo = limt-xx, (/(*)) = lim t _ >00 (/(x,t)) x « p(t = 1800). This 
approximation only breaks down extremely close to the phase 
boundary between the red and green regions and the differ- 
ence is not visible at the resolution of the plot. The dotted 
line shows the mean field transition. The dashed line shows 
the approximate position of the directed percolation phase 
transition for a linear expansion, as shown in Fig. [JJ 



computed using Laplace transform techniques (see e.g., 
Ref. 03, H|) and is given by: 



S r (r) 



2tt 7r * — ' n 

n—l 



00 



e -AnV/4 sin 



n<j) 



(49) 



Note that the results for the survival probability in 
Eq. (|4U)) and Eq. (|4T|) coincide when A is small, and the 
finite range of the sector size becomes irrelevant. 



D. Directed Percolation With Inflation 

As illustrated in Fig. [TJ linear range expansions within 
the Domany-Kinzel (DK) model on a hexagonal lattice 
for /if, = lead to a phase transition driven by com- 
petition between the selective advantage s of the green 
cells and deleterious forward mutations \i / , from green to 
red. In population genetics, a two allele model with ir- 
reversible mutations is an important limiting case in the 
theory of quasi-species [39144 1| . In this theory, the phase 
transition corresponds to an "error threshold" at which 
a well adapted population's fitness distribution near a 
single fitness maximum becomes delocalized due to the 
large phase space available for deleterious mutations (see 
Appendix |A"1 for more details). In this section, we will 
compare the well studied DK model results on a hexago- 
nal lattice [13, [24|, [2?1 [HI with our radial range expansion 
model results on a disordered Bennett lattice. 
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When selection is strong, the green strain is able to sur- 
vive even when p,f > 0. However, for large enough pj, 
we expect a "mutational meltdown" such that the advan- 
tageous green strain eventually goes extinct [HI, [2(| . It 
is known that in the biologically relevant region of rea- 
sonably small, nonzero deleterious mutation rates (e.g. 
< fif < 0.5), the model exhibits a phase transition in 
the directed percolation universality class. This transi- 
tion is in the one-dimensional compact directed perco- 
lation (CDP) or "voter model" universality class when 
M/ = Mb = 0- The phase diagram as a function of param- 
eters analogous to our parameters s and pf for /i& = 
has been studied extensively (see (2(| HI] and references 
therein) and is illustrated in Fig. [TJ Note that the CDP 
transition is connected to an entire line of directed per- 
colation phase transitions. 

Near the phase transition line, various dynamical 
quantities as power laws in time, similar to the 1/t be- 
havior of the density of surviving green cells with random 
initial conditions found in mean field theory (Eq. ([AT])) 
[20T [H]. In addition, for a fixed mutation rate pf > 0, 
one can show that there are four critical exponents that 
characterize the scaling of the system. Two exponents, 
v±_ and , describe the diverging width and length of 
the growing green clusters parallel and perpendicular to 
the time-like direction. If the directed percolation tran- 
sition occurs at a critical selective advantage s c (/i/), and 
5 = s — s c (0 < 5 <C 1) measures the distance from the 
critical point in the active phase, then a typical width 
and length of the green sectors diverge as £j_ ~ 8~ u± and 
£|l ~ 5~ v n , respectively. 

Two additional exponents, j3 and /?', describe order pa- 
rameters relevant for two distinct initial conditions: (1) 
the average density of active sites (green cells) at long 
times, poc = (f(t — > 00)) starting from /(0) = 1; and (2) 
the long-time sector survival probability Sqo = S ss (t — !> 
00 ) of a single green cell (more generally, any small green 
sector) seeded in an all red homeland. Here S SB (t) is the 
probability that a single sector has survived out to time 
t. We then have poc ~ 5@ and ~ for small dis- 
tances 5 > into the active phase. In the biological con- 
text, the exponent /3 describes the mutational meltdown 
in a population of green cells during a range expansion, 
and P' characterizes the survival probability of a rare ad- 
vantageous mutant that spreads through the population 
while experiencing deleterious back mutations to a less 
fit strain. See Fig.[TTJfor Poo{pf, s) for radial expansions. 
Although we don't expect a sharp phase transition (see 
below), the general characteristics mimic the linear ex- 
pansion diagrammed in Fig. [TJ 

An important feature of the directed percolation (DP) 
dynamics described by Eqs. (0 fTUj) is the so-called ra- 
pidity reversal symmetry that implies /3 = j3'. Rapid- 
ity reversal is a special kind of time reversal that can 
be seen explicitly in the field theoretic formulation of 
the Langevin equation, which is known to be equiva- 
lent to a Rcggcon field theory [2(| [28| . Rapidity rever- 
sal symmetry is only valid asymptotically as S — > + , 




FIG. 12: Linear and radial range expansions with pb = and 
Pf = 0.1 for various s near the directed percolation phase 
transition. Green organisms have a selective advantage s, but 
are subject to forward mutations to a less fit red strain. In 
the radial expansions (a) Rq = 50 and the system evolved 
for ~ 250 generations. The dark blue circle indicates the 
crossover time t* = 100. In (i>) a linear range expansion of 
300 cells evolves for 400 generations. 
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FIG. 13: Data collapse for a radial range expansion with 
Ra = 2048 for the average number of green cells over about 
600 generations. We test many values for r = s — s c = 
±0.01,0.02,... for p — 0.1. The collapses are consistent 
with the directed percolation values a m 0.159464(6) and 
~ 1.733847(6). The critical selection parameter was ap- 
proximately s c ~ 0.55(1). 



where it implies that p{t) = (f(t)) oc S ss (t) near the 
directed percolation phase transition for large values 
of t. Indeed, both quantities scale for long times like 
p(t) - S ss (t) - t~ a , where a = fi/v\\ = /v\\. The 
asymptotic relation p(t) cx S ss (t) is also valid for systems 
with time-independent, small system sizes since the size 
is invariant under time reversal (see Ref. (23j ] for an exam- 
ple). However, unlike a typical finite size effect, inflation 
is inherently asymmetric in time, and our simulations 
will show that inflation breaks the rapidity reversal sym- 
metry, leading to new physics beyond conventional finite 
size scaling. 
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Consider the DP transition in the radial setting. The 
plot in Fig. [TTJ of the averaged green cell density after 
1800 generations (p(t = 1800) « poo) for an all green 
homeland with Rq = 200 shows that we have a crossover 
between an active and an inactive phase. Despite the ex- 
istence of these two phases, interpreting the critical point 
for radial expansions as a DP phase transition requires 
care! Inflation inhibits the boundaries of green sectors 
that have formed during time t* from interacting diffu- 
sively at longer times t 3> t* ; hence a large enough green 
sector could, with small probability, survive to the infla- 
tionary regime even in the inactive phase, after which it 
will persist for a very long time. Of course, increasing the 
forward mutation rate pf still creates a crossover in the 
radial expansion dynamics, since any green sector will 
eventually become contaminated with red mutations. 

The crossover in a radial range expansion is slightly 
shifted relative to the linear expansion transition line 
(dashed line in Fig. ITT)) . The precise value of the transi- 
tion is not universal and varies with the lattice model 
details. The direction of the shift in our case might 
be explained by the larger effective coordination number 
z r > 2 in the disordered Bennett lattice relative to the 
hexagonal lattice: A larger z should move the transition 
closer to the mean field line (dotted line in Fig. ITT)) as 
each daughter in a new generation can evolve from more 
potential parents. 

For early times t < t', inflation should only have a 
minor effect on the coarsening dynamics at the transition. 
As a check, Fig. [T3] shows a collapse of the data for the 
average fraction pit) = (f(t)) of green cells for a variety 
of selective advantages starting with an all green initial 
population with Rq = 2048. Our simulations run for 
about 600 generations for different values of s at a fixed 
pj = 0.1. We successively collapsed the data using the 
remarkably precise estimates for the directed percolation 
exponents a w 0.159464(6) and v\\ w 1.733847(6) quoted 
in Ref. [2Cj. As expected from Fig. [TTJ the position of the 
apparent critical point for times t <C t* is shifted slightly 
relative to the linear case (s c w 0.55 versus s c w 0.57 for 
(if = 0.1). Fig. 1 131 confirms that the early time dynamics 
of our model is indeed accurately governed by directed 
percolation. 

At late times t S> t* , the population fronts become lo- 
cally flat due to the increasing population radius. Thus, 
one might naively expect that the f > f dynamics at 
the transition is also governed by critical DP dynamics. 
However, inflation causally disconnects portions of the 
population after time t*. Thus, the J > f dynamics 
consists of a set of noninteracting DP processes. The 
fluctuations around time t* set the initial conditions for 
each process. These initial conditions can be quite differ- 
ent (ranging from all green cells to a single green cell) and 
lead to different DP dynamics with different power laws 
j20t [28j . Since we average over the behaviors of all the 
sectors in a radial range expansion, the t 3> t* regime 
does not correspond to the usual DP dynamics with a 
single initial condition. 



One case in which the long-time behaviors of the linear 
and radial range expansions are similar is in the limit 
Pb <C s for an all green initial condition, in which the 
expansions both achieve an active steady state. The dy- 
namics is deep in the "active" phase - mutations generate 
small red sectors which die out with high probability. If 
these sectors are sufficiently dilute, we can neglect their 
interactions. This phase is illustrated on the right panel 
in Fig. I12r a). Selective advantage bias tends to squeeze 
out red sectors, which we model by first replacing w by 
— w in Eq. pip to study the probability distribution of red 
sector widths. However, unlike the analysis of green sec- 
tor widths in Sec. IIII CI the "homeland" radius Rq must 
now be replaced by R(t r ) = Rq + vt r > Rq, where t r is 
the time at which a red sector is generated by mutations. 

To obtain the long-time behavior of the density p^ = 
p(t — > 00) = (f(t — > 00)) of green sites, we examine the 
survival rate of red sectors formed at very large times t 3> 
t* . In this case, the fate of the red sector is only weakly 
influenced by inflation. The red sector size probability 
distribution then obeys 



_ dp <9 2 p 

d tPs (<j>,t)=w-+A w , 



(50) 



with the 



istant coefficients, for a fixed time t' > t*. 



7 R(t) 
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,(R(tW 



(51) 



As usual, we impose an absorbing boundary condition 
p(4> = 0, i) — on Eq. (|50l) . If we change variables to 
make Eq. ([50]) a conventional diffusion equation, then the 
boundary condition now advances on the domain where 
diffusion and selection bias take place. In such a case, 
the absorbing state (extinction of a red sector) is always 
reached [26[. Finally, note that we arc no longer strictly 
in the small s regime and hence have no simple expression 
for 7 in Eq. ([51]). However, we expect it to be close to 
unity (see Appx. [D]) . Since we are only interested in an 
approximation of the scaling of (f(t — > 00)) with pf and 
s, we set a r = 7 = 1 in the following for simplicity. 

The average excursion area (red sector size) (A rc d) of 
the random walk described by Eq. ([BTJ[) before it reaches 
the absorbing state, to first order in the initial sector size 
of one cell, (f>o ~ a r /R(t), is given by a known result from 
diffusion theory (see e.g., Ref. [42j]): 



(Ared) — 7=\2 



AR(t)(j)Q 



1 



(52) 



Next, consider a population of green cells (in the active 
state) evolving from time t = t to time t + T (with t ^ 
T,t*). The number of red sectors iV rG d that are seeded 
in the population during this time will be approximately 
A rc( j w 2irp,fvTR(t)/ar- The total area covered by the 
population from time t to time t + T is A tot w 2ixR(t)T. 
Upon dividing the average area covered by the green cells 
by Atot, we hud the steady-state density of the green 
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FIG. 14: We plot the average green cell density p(t) = (/(£)} 
at long times t versus a control parameter Pf/s 2 for an ini- 
tially all green population for radial and linear expansions. 
For radial expansions, f(t) is averaged over 10 3 runs at 
t = 2 x 10 3 generations. In the linear case, we averaged 
f(t) over 2 x 10 3 runs, at time t = 2 x 10 4 , for a system size 
of 10 3 . Away from the transition region (pf/s 2 < 0.3), the 
time t is long enough so that the density p(t) for both ex- 
pansions is approximately equal to the steady-state density 
Poo = {f(t — > oo)}, with an error smaller than the point size. 
The black line shows the predicted scaling of the steady-state 
density for small pf/s 2 . We varied s between and 1 and pf 
between and 0.5. 
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Note the crucial difference between Eq. (|53|) and the 
steady state value 1 — fXf/s in the mean field (well mixed) 
solution in Eq. (|A1[) : For fixed values of pif, we require 
larger values of s rj J~pTf > pf to get a comparable re- 
duction in p for the circular range expansion. 

Eq. (|5"3"]) is derived in the absence of inflation, so the 
same expression for p holds for linear range expansions, 
similar to the result obtained for a related model in 
Ref. [18|. This is consistent with the similarity between 
the active phase pictures for radial and linear range ex- 
pansions in Fig. WMa) and (6), respectively. We now 
check Eq. ([5"3"]) against our data in the phase diagram of 
Fig. [TT] for both linear and radial expansions. Fig. [HI 
shows collapse of the data for p,f/s 2 <?C 1, i.e., deep in 
the active phase. As argued in Ref. [Hj| . p is driven be- 
low 1 — Pf/s 2 by collisions between red sectors, which 
merge to shield green sectors from invading the red re- 
gions, increasing the total fraction of red cells. Red sector 
collisions are suppressed by inflation, so p stays closer to 
the 1 — ^f/s 2 line for radial expansions as opposed to 
linear ones, as is evident in Fig. [TU Note that the slope 
of the line at small Pf/s 2 in Fig. [14] is consistent with 
our approximation a r = 7 = 1. 

As we keep increasing Pf/s 2 , the steady-state density 
Poo eventually goes to zero. In the linear case, there 
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FIG. 15: Data collapses for radial range expansion simulations 
with many values of t* — Ro/v (with pj =0.1 and s = 0.547). 
In (a) we plot the scaled green cell density p(t) — (f(t)) with 
an all green initial condition for Rq = 8, 16, 32, . . . , 2048. The 
regimes t <C t* and ( > (* can be described by power laws. 
The t~ a decay for t < t* with a « 0.16 is expected from 
directed percolation theory [2(j, but the exponent a ~ 0.26 
that arises for long times appears to be new. In (b) we plot 
the scaled sector survival probability S ss (t) of a single green 
cell inoculated in an all red homeland Ro — 16, 32, . . . , 2048. 
The dashed line indicates the expected directed percolation 
result: a decay with slope —a ~ —0.16. Note, however, that 
the survival probability levels off to a finite value for J > (* 
due to inflation. 



is a sharp transition at around p-f/s 2 « 0.3, consistent 
with the results of Ref. [IH . This extinction of the green 
organisms is triggered, of course, by crossing the directed 
percolation phase transition line [201 ] . A sharp transition 
is less evident in Fig.[14]fbr radial expansions: the density 
p seems to decay to zero more slowly as we approach the 
crossover. We shall now study the dynamics of the radial 
range expansion near this crossover. 

To explore directed percolation dynamics when infla- 
tion has set in, i.e. for t > t* = Rq/v 7 we use a dy- 
namic scaling hypothesis for the average green fraction 
p(i) = (f(t)) and for S ss (t), where S ss (t) is the survival 
probability of a single green cell inoculated at the edge of 
a green homeland. We let 5 = s — s c (pf) be the distance 
from the linear inoculation (i?o — > 00) critical point. We 
approximate the critical point s c (/x/) by looking at the 
short time (t <C t*) dynamics of a radial range expan- 
sion with a large homeland radius Rq. For example, we 
found that s c (p f = 0.1) w 0.55 (see Fig. [□]). Note that 
we cannot use the hexagonal lattice for linear expansions 
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to find the critical point because we expect the point to 
shift when we use the disordered Bennett lattice. 

Inspired by finite size scaling ideas near critical points, 
we introduce a scaling transformation S — > bS and expect 
the scaling relations 

p(t) = b- p p [&-"« t, bS, b v t*\ , (54) 

s ss (t) = b-P's ss [b- v n,bs,b- v n*] , (55) 

where we have suppressed metric factors [20| . and p and 
5 SS are scaling functions. We now set b~ u wt* = 1 and find 
that our data for various t* (corresponding, say, to differ- 
ent homeland radii Rq — vt*) at the critical point (<5 = 0) 
should collapse upon plotting p(t)(t*) a and S ss (t)(t*) a 
versus t/t*. These scaling forms are tested in Fig. [T5] 
The data collapse reasonably well, and we indeed get 
a crossover at t/t* « 1 in both cases. The inflation 
crossover time t* enters just as a finite size effect would 
in conventional critical phenomena. 

Both of the collapses are consistent with a single value 
a given by directed percolation theory [20j . This suggests 
that we retain the rapidity reversal symmetry of directed 
percolation in a linear geometry at early times. How- 
ever, the rapidity reversal prediction |20| p(t) cx S ss (t) 
is violated in the inflationary regime t > t* . Inflation 
prevents green sectors from dying out for a single seed 
initial condition (Fig. ITSb). and the survival probabil- 
ity S ss (t) approaches a nonzero constant at long times. 
Conversely, inflation evidently enhances the loss of the 
total fraction of green cells for the green homeland initial 
condition, and the density p(t) decays to zero as t — > oo 
(Fig. In a conventional finite size effect, S ss (t) and 

p{t) are consistent with rapidity reversal symmetry, and 
both decay exponentially to zero at long times [231 ] . 

The levelling off of the survival probability S ss (t) in 
Fig. IT5T 6) for t ^s> t* is consistent with the analysis of 
single sectors without mutations in Sec. MI CI The sur- 
vival probability of a single sector approaches a constant 
at long times (see Eq. |4"5]) even at s = 0. Unlike a sector 
in a linear expansion, the diffusive motion of the sector 
boundaries in a radial range expansion is suppressed in 
the inflationary regime. The boundaries are unlikely to 
collide in this regime, and a sector that makes it to the 
inflationary regime can survive even as t — > oo. In a 
linear expansion, S ss (t) decays to zero at the directed 
percolation transition as t~ a for large times t. 

The more rapid decay of the average fraction of green 
cells p(t) in Fig. [IBTa) for t > t* is more subtle: p(t) 
seems to decay approximately as t~ - 26 . We can under- 
stand this trend qualitatively by appealing to the mean 
field solution for p{t) (see Eq. (jAllO which decays at the 
much faster rate t~~ x at criticality. Since inflation sup- 
presses both diffusion and noise in our Langevin equation 
for /(x,i) (Eq. (UU)), it is plausible that p(t) = (f(t)) de- 
cays faster than for the linear expansion where diffusion 
and noise are more important. However, it is not clear 
that the larger decay exponent in the inflationary regime 
is universal, and we do not have an explanation for its 
particular value of —0.26. 



Another important dynamical quantity is the pair con- 
nectedness correlation function [28j |. i.e., the probability 
T(x,t;xo,to) of finding a directed, causal chain of en- 
tirely green cells between a green cell at point (xo,io) 
and any point [x, t) for t > to. For radial Domany-Kinzcl 
expansions, T(x, t\ xq, to) is the probability that a green 
cell at position xq = R(to)<j>o along the population cir- 
cumference at time to has a green cell descendant at some 
later time t at position x — R(t)<j>. Thus, T(x, t; Xo, to) 
characterizes the spatial correlations of genetic lineages 
in the green cells. The pair connectedness function does 
not keep track of the red cell lineages, including those red 
cells that evolved from green cell ancestors. A knowl- 
edge of T(x,t;xo,to) might be useful for genetic infer- 
ence, where spatial distributions of genetic variants at 
time t are used to infer ancestral genetic distributions 
(see Ref. [l(| for more information). 

We can probe T(x, t; xo, to) by employing a single green 
cell initial condition at the origin and monitoring the 
resultant green sector formed by descendants of the single 
green ancestor. We compute the total number of green 
(i.e. active) cells N a (t) at time t and the mean square 
spread of green clusters that survive to time t: X%(t) = 
(|Ax| 2 ), where 



(|Axp 



J(Ax) 2 T (Ax,£;0,0) d(Ax) 
/T(Ax,i;0,0)d(Ax) ' 



(56) 



and we average over all green cells located at displace- 
ments Ax = R(t)A<fi away from the initial green cell. 
The quantity N a (t) is related to T(x,i;0,0) via [H 



N a (t) = const, x J dx T(x, t; 0, 0) 

and obeys a scaling relation (see Ref. (20l. [28^) 

N a (t) = b^-' 3 -' 3 ' 'N a (b- u H,bS,b~ u H*). 
Similarly, we have 

const. 



xi(t) 



N a (t) 

and the scaling relation 



dx|x| 2 T(x,t;0,0) 



x 2 (t) = b 2 ^R 2 a (b-^t,bS,b-»n*), 



(57) 



(58) 



(59) 



(60) 



where we have again suppressed metric factors [28[. 

Upon setting b~ u wt* = 1, we now find data 
collapse when we plot X%(t)(t* )~ 2v J-/ v \\ and 
N a {t){t*)^+P'- V ^/^ versus t/t* at criticality (S = 0) 
(see Fig. [T6|) . In our simulations, X 2 (t) can be computed 
by finding the average squared angular spread (<f> 2 ) and 
multiplying by R 2 (t) = (Rq + vt) 2 . We find excellent 
data collapses for both quantities. Notice that for times 
t <C t* , we get the expected power- law behavior given 
by the directed percolation universality class (28|. 

For t ^> t* , we have argued that a green sector evolved 
from a single green cell will expand dctcrministically. The 
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FIG. 16: Data collapse for radial expansion simulations with a 
single green cell in an otherwise red initial homeland perime- 
ter (with (if = 0.1, fit = 0, and s = 0.547) on log-log 
plots, (a) The rescaled mean squared cluster spread X%(t). 
The dashed line for t <C t* indicates the directed percola- 
tion power law scaling with exponent: 2v±/vu fa 1.3. The 
dashed line in the inflation-dominated t 3> t* regime agrees 
with a Ra{t) ~ t 2 behavior discussed in the text, (b) Plot 
of the rescaled average number of green cells N a (t). The ex- 
pected directed percolation scaling for t <C t* is shown with 
a dashed line, with (u± — 2fi)/v\\ ~ 0.31. For t 3> t*, we 
find N a (t) ~ f 1_Q ~ i ' 84 , consistent with arguments incor- 
porating inflation given in the text. The simulation consisted 
of 9.6 x 10 4 independent runs for each homeland with radii 
Ro = 10, 30, 50, 100, . . . , 2000, evolved to a final population of 
size R(t) = 2800. 



angular width of the sector will increase somewhat due 
to the selective advantage, but inflation will be the dom- 
inant effect at long times. Hence, the main contribution 
to X 2 a for t > t* arises from large green sectors evolving 
with fixed angular widths. We then expect X^(t) ~ t 2 
at long times, consistent with the long-time behavior in 
Fig. [ToTa). The main contribution to N a (t) also arises 
from this deterministic expansion epoch. However, we 
expect the density of green cells inside each sector to de- 
cay according to t~ a , since the interior of a green sector 
should be describable as a critical directed percolation 
process. Thus, we expect that N a (t) ~ t 1 ~ a for t > t*. 
This expectation also matches the simulation results (see 
Fig. [T6J&)). By contrast, directed percolation in a fi- 
nite geometry is unable to sustain large cluster growth. 
Thus, at long times, N a (t) and X 2 (t) have different scal- 
ing functions that rapidly decay to zero [23j . 

The scaling result for the mean square spread X 2 of 



FIG. 17: An inoculation of two bacterial (Escherichia coli) 
strains on a Petri dish, labelled with blue and yellow consti- 
tutively expressed fluorescent proteins. The strains are oth- 
erwise identical and, consequently, represent the s = fj,f = 
fib — case in our analysis. The initial inoculation has the 
shape of an partial ring to allow for nutrients to penetrate the 
ring center and allow for the bacteria to grow inward, creat- 
ing a "deflating" radial range expansion. A regular "inflat- 
ing" radial range expansion moves away from the ring. This 
single experiment illustrates both kinds of radial expansion 
discussed in this paper! 



surviving green clusters suggests a simple prediction for 
directed percolation with arbitrary power law inflation 
R(t) = R [l + {t/t*) e }. Namely, for 9 > u±/u\\ w 0.632, 
the inflation should be able to overtake the DP critical 
cluster growth, so that X 2 ~ t 26 and 7V a ~ t e ~ a at long 
times t ^> t* . X 2 and N a should still have their regular 
DP power laws for t <C t* . However, the shape and 
position of the crossover region at t ~ t* will likely change 
with O. Conversely, when < v±/vi\, the cluster growth 
is always faster than the inflation, and we expect the 
inflation to be irrelevant so that X 2 ~ t 2 ^/"!! and N a ~ 
t (u ± -2/3)/u t at long times for all q < v± / v ^ Of course, 

a rigorous check of this simple prediction is necessary. 
It would also be interesting to study the borderline case 



IV. RADIAL RANGE EXPANSIONS WITH 
DEFLATION 

Another type of radial expansion starts with a popu- 
lation around the circumference of a circle of radius Ro. 
The population then grows inward toward the center of 
the circle, colonizing the interior. An example of such a 
"deflating" range expansion is the settling of an island 
by a new species that arrives simultaneously around the 
shoreline and then grows inland. The population front 
at time t in this case is a circle of radius R(t) = Ro — vt. 
Note that t* = Ro/v is now a sharp temporal cutoff (in- 
stead of a crossover), at which the population front size 
vanishes and the invading species completely takes over 
the circular island. 
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a = 0.01 




FIG. 18: (a) A radial range expansion with a deflating pop- 
ulation front, for s = fif = fj,b = 0, starting from an initial 
population around a circular rim of radius Ro = 300. We use 
the same Bennett model lattice described in Sec. [Til while run- 
ning the update scheme in reverse, (b) Deflationary range ex- 
pansions for a single green sector with an initial size <f>o = 0.1 
radians for a large (s = 0.4) and very small (s = 0.01) selec- 
tive advantage. In both cases, deflation forces the green sector 
boundaries to collide at the center (when t = t* — Ro/v). 



A deflating range expansion can be realized in experi- 
ment by inoculating a ring of organisms on a Petri dish, 
as shown in Fig. [TTJ As long as there are enough nutri- 
ents on the interior of the ring, the organisms will grow 
inward toward the ring center. Fig. [17] shows the ge- 
netic demixing occurring between two neutral bacterial 
strains labelled with constitutivcly expressed blue and 
yellow fluorescent proteins. The bacterial sectors exhibit 
a counterclockwise "twisting" that is analyzed in detail in 
Ref. [l[. This uniform tangential motion of the bacterial 
sector boundaries does not affect the genetic demixing or 
quantities such as the heterozygosity that only depend on 
the separation between two individuals along a popula- 
tion front. We will not model the twisting in simulations, 
but note that it can be easily implemented by shifting all 
of the cells at a population front counterclockwise at each 
time step. 

The simulated evolution of two neutral strains in a de- 
flating range expansion is shown in Fig. H8f o). A well 
mixed initial condition is used with equal proportions 
of green and red cells. The initial population is on a 
circular boundary with radius Rq = 300. The deflation- 
ary dynamics, similar to the inflationary dynamics con- 
sidered in Sec. IIII A[ is related to a linear range expan- 
sion via a conformal time coordinate t c (t) = t/(l — t/t*). 
The heterozygosity is again given by H(<f>, t) = Ho erf |£|, 
with £ = Rq4>/ y/8D r t c (t). However, unlike the inflation- 
ary case, H((j>, t) now decays to zero and becomes flat as 
t — >• t*\ Thus, the deflating range expansion evolves to 
a population with just one genetic sector as t — > t* , as 
shown in Fig. ITST a). 

The dynamics of a single green sector in a deflating 
radial range expansion (shown in Fig. [PSt^ )) also differs 
from the inflating case. Using the diffusion equation tech- 
niques of Sec. IIII CI and App. O we can calculate the 




FIG. 19: A comparison of the long-time survival probabil- 
ities of green sectors for small selective advantages w oc s 
in range expansions with radial inflation (S r {t — > oo), solid 
lines), radial deflation (Sd(t — > t*), dashed lines), and for 
linear range expansions (Si(t — > oo), dotted lines). The 
probabilities are calculated using Eq. (|45p . Eq. (|6ip . and a 
similar result for a linear range expansion. In the forward 
and reverse radial expansions, the initial population radius is 
Ro = 100 a r , where the effective lattice constant a r — 1 for 
convenience. The linear expansion assumes a total population 
size L = 2tvRo = 200-™,-. We used two different initial green 
sector sizes, 0o = 0.01, 0.002. 



deflationary survival probability Sd(t) of a green sector 
with a small selective advantage w and a small initial size 
00 < 27T. The result as t -> t* is 



S d (t->t*)«l 



dx, 



(61) 
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Note that the convergence of sector boundaries insures 
that Sd(t — > t*) vanishes in the absence of selection 
(w = 0). This behavior differs markedly from the in- 
flationary case in which S r (t — > oo) has a limiting value 
given by Eq. (|47[) with t = 1. However, when t — > t* 
in the deflationary case, the population size at the fron- 
tier goes to zero, and the finite size of the cells becomes 
important. Hence, we do not expect the continuum dif- 
fusion model to hold exactly as t — > t*; contrary to the 
continuum result (Eq. (|6ip ). a small residual survival 
probability is possible in this limit, especially for larger 
initial sector sizes 4>o- Thus, Eq. (|5Tj) and Eq. (|52"j) are 
good approximations to the limiting survival probabil- 
ity only when the fate of the green sector is determined 
sufficiently far from the sector convergence time t = t* . 

Comparing Eq. (|6"2"|) and Eq. (f4"gj) reveals that defla- 
tion suppresses the effects of selection relative to infla- 
tion, an effect which arises because deflation forces the 
green sector boundaries closer together, as illustrated in 
Fig. I18f &) . The survival probabilities for inflationary and 
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deflationary radial range expansions can also be com- 
pared to the limiting survival probability Si(t — > oo) of 
a green sector in a linear range expansion with selective 
advantage w and initial size i?o</>o in a population of lin- 
ear size L = 2ttRo. The results for two different sector 
sizes <po are shown in Fig. 1191 As expected, the survival 
probability of a green sector in a linear range expansion is 
bracketed by the probabilities for inflating and deflating 
radial expansions. 



V. CONCLUSIONS 

This paper has explored differences between the pop- 
ulation genetics of radial and linear range expansions in 
the presence of genetic drift, mutation, and selection. We 
find that the inflation embodied in a radial expansion de- 
creases the effect of genetic drift and biases the stochastic 
dynamics in a deterministic direction. The proliferation 
of extra sites at the outer edge of the radial expansion 
greatly reduces sector interactions after time t* = Rq/v. 
For times i 3> t* , the population evolves in isolated angu- 
lar segments, with boundaries perturbed by genetic drift. 
Near the directed percolation (DP) phase transition, we 
found that t* acts as a finite crossover time, after which 
the radial expansion no longer experiences the critical 
dynamics. Finite size scaling near the DP transition due 
to inflation is described by a different scaling function 
than that governing conventional finite size effects. 

We presented and tested analytical approximations to 
the dynamics of single sectors in various limits, includ- 
ing selective advantages weak compared to the genetic 
drift and inflation. It would be interesting to develop a 
more precise asymptotic analysis of the radial survival 
probability S r ((j>o,T) of a sector with an arbitrary initial 
angular width </> and r = t/(t + t*). We also calculated 
the survival probability in a deflationary radial range ex- 
pansion in which the sector boundaries are pushed to- 
gether by a collapsing population front. We found that 
inflation enhances and deflation diminishes the survival 
probability relative to a linear range expansion with the 
same initial population size and selective advantage. 

It would also be interesting to study the inflationary 
effects in more detail via a field theoretic derivation of the 
scaling hypothesis [2(| [28| . In the linear expansion, the 
relevant field theory is Reggeon field theory. In the radial 
case, we would have to introduce the inflationary metric 
ds 2 = v 2 dt 2 + R(t) 2 d<f> 2 (with front propagation speed 
v = dR/dt) for measuring distances between points in 
the range expansion rather than the usual Euclidean one 
d.s 2 = v 2 dt 2 + dx 2 to take into account the expanding 
system size. So, we expect that the field theoretic de- 
scription of radial expansions would involve studying the 
Reggeon field theory in this "curved" inflationary space. 

It would also be interesting to extend our results by 
allowing for undulations in the radial population front. 
Such undulations produce small changes in the power 
laws predicted by simple diffusive models of sector growth 



and wandering [18| . One way to allow undulations in 
simulations is to generalize the Eden model proposed 
in [l9[ to an amorphous lattice and evolve initially cir- 
cular homelands of individuals which now divide freely 
at the population periphery. Without the constraint 
that daughter cells always spawn as close as possible to 
the center of the expansion, the population front should 
roughen away from a nearly circular shape. In the fu- 
ture, we hope to study an extension of our model to 
three dimensions, thus allowing analysis of a greater va- 
riety of range expansions, such as tumor growth in can- 
cer. Spherical cell colonies, for example, can be modeled 
using the three-dimensional Bennett model (2{J. We ex- 
pect that the effects of inflation are more pronounced 
in three dimensions since sector surface areas now in- 
crease quadratically in time due to inflation, while sector 
coarsening due to genetic drift is much slower at a two 
dimensional frontier II Oil - 
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Appendix A: Quasi-species and Mutation-Selection 
Balance 

An important biological motivation for studying irre- 
versible, deleterious mutations is the molecular quasis- 
pecies theory developed by Eigen (39[. In this theory 
(reviewed in [4(| El| ) , each individual in an asexually re- 
producing population has a molecular sequence (i.e., a 
DNA or a protein sequence). During each reproduction 
step, the sequence gets replicated and acquires errors at 
some rate [i. The errors change the sequence's fitness. 
If the error rate is small, the population should become 
well adapted and its sequences will be clustered around a 
peak in a fitness landscape. The individuals at the peak 
have a unique "master sequence" corresponding to this 
highest possible fitness in the population, as illustrated in 
Fig. 1201 Each mutation away from the master sequence 
is necessarily deleterious. For a variety of fitness peaks 
with this general structure, there is a critical mutation 
rate /x c above which the population fitness distribution 
becomes delocalized, and the fraction of individuals that 
preserve the master sequence goes to zero at long times. 
The critical mutation rate [i c is called an "error thresh- 
old" 0. 

The irreversible mutation case (fib = 0) we consider 
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• cell with master seq. 

• cell without master seq. 



d(f(t))/dt=(s-(if)(f(t))-s(f(t)) 



2 . 




sequence space 



sequence space 



FIG. 20: A schematic of the quasi-species model in which 
each cell in a population has a sequence distributed in a fit- 
ness landscape with a single sharp peak (in green) surrounded 
by a flat landscape (in red). As cells divide, each individ- 
ual reproduces with a mutational error rate \i that moves it 
around in sequence space (black arrows). Below a critical 
"error threshold" (fi < /i c ), a finite fraction of the population 
has the master sequence (i.e., is a green cell). Cells deviating 
from this ideal sequence are red. When (i > (i a , the popula- 
tion fitness distribution becomes delocalized, and the fitness 
peak is depopulated. 



in Sec. IIII Dl corresponds to a particularly simple fitness 
landscape in which the fitness peak is infinitely sharp and 
each mutation away from the master sequence will have 
an identical deleterious effect with a selection coefficient 
s (see Fig. [2U] for an illustration and Ref. for more 
details) . The green cells are the cells with the master 
sequence. Cells without the master sequence are red and 
have the same fitness because the landscape is flat away 
from the peak. If the master sequence is very long, then 
the probability that a red cell recovers the master se- 
quence as it wanders aimlessly through sequence space 
is vanishingly small. We can then treat the error rate (i 
as an irreversible mutation rate (i = (if which mutates 
green cells to red cells. The red cells also acquire errors 
with rate (i, but these errors will just move them around 
the flat fitness landscape. 



The error threshold transition in the quasi-species 
model for the simple fitness landscape in Fig. [20] is iden- 
tical to the directed percolation transition discussed in 
this paper. In population genetics, these transitions are 
typically studied in well mixed populations, which cor- 
responds to a mean field solution in the directed perco- 
lation language. The mean field solution for the den- 
sity of green cells f(t) (fraction of cells with a master- 
sequence) can be found from Eq. © by assuming a ho- 
mogeneous (space independent) density /(x, t) = f(t) 
and negligible noise. For an all green initial popula- 
tion (f(t = 0) = 1), the fraction of green cells (aver- 
aged over many experiments) is given by the solution of 
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Thus, (f{t)) approaches Oor 1 — fif/s exponentially with 
time for (if < s and (if > s, respectively. The er- 
ror threshold is given by the critical selective advan- 
tage s c = (if, and we evidently have a phase transi- 
tion controlled by mutation-selection balance. For more 
details about the relationship betwee n q uasi-species and 
mutation-selection balance, see Ref. [4J]. Note that this 
mean field analysis assumes a very large population size 
with negligible genetic drift - the sharp phase transition 
is destroyed by number fluctuations in the well mixed 
"zero-dimensional" case. 

Finite spatial diffusion and noise due to genetic drift 
drastically change the position of the phase transition 
and the broad features of the dynamics of /(x, t). Upon 
defining the rates (i / and s in units of the inverse cell divi- 
sion time (equivalent to setting the generation time r g = 
1), we construct the phase diagram (Fig. [TTj) by study- 
ing the limiting fraction of green cells lim t _ >00 (/(x, t)) x 
with the "all green" initial condition /(x, 0) = 1, where 
(. . .) x is both an ensemble average and an average over 
all cell positions in the growing population at some time 
t. The Monte Carlo simulations used to calculate this 
quantity are described in Sec. [TXJ A much stronger crit- 
ical selective advantage (s c ~ -JW}) than that expected 
from mean field theory (s c = (if) is required to main- 
tain the green strain in the population. These and other 
corrections to the mean field behavior have been exten- 
sively studied with both simulations and field theoretic 
methods [11 El. 



Appendix B: Coarse-Graining and Langevin 
Description 

The simulations in this paper focus on voter-type mod- 
els, with one individual per lattice site. However, a 
coarse-graining procedure leads to a description in terms 
of continuous variables, similar to that reviewed for step- 
ping stone models in Ref. We coarse-grain voter- 
type models by considering a set of 51 > 1 cells, and 
rcscaling our lattice constant so that these fl cells are 
all contained within a spatial volume a d . Following, e.g., 
Vazquez and Lopez [44| , a coarse-grained population den- 
sity of green cells is given by 



l + o-,-(t) 



n 2 



(Bl) 



where the sum is over cells in a neighborhood Sq(x) of 
the fl cells centered on x, and we use Ising spin variables 
Ui{t) = ±1 to specify the occupancy of site i. Individual 
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cells in this neighborhood change states according to the 
rates given by Eq. Q. We now randomly choose one of 
the fl spins to flip in each neighborhood. This means 
that each spin flip at x' changes our field /(x, t) by a 
displacement field ±r(x,x') given by 



±r(x,x') = ±fi-V<5(x - x') 



(B2) 



where the sign depends on the state of the spin at x'. In 
addition, we now use the average value of the spins in 
neighboring cells to evaluate the summations over near- 
est neighbors in the flipping rates due to selection and 
genetic drift (see Eq. © and Eq. (|7|)). The resulting 
dynamics resembles the stepping stone model of popula- 
tion genetics [HI IH| , with populations of size ft on com- 
pact "islands" where mutation and selection act, with 
exchange of individuals between neighboring islands. 

To study the dynamics of the the coarse-grained 
population density /(x, t), we introduce rates 
w(/(x); ±r(x, x')) that describe a transition from a 
starting field configuration / to a field fir. The new 
coarse-grained transition rates (in units of t^ 1 ) that 
define an underlying master equation for the field / are 



w[/;+r] = w[/(x)->/(x)+r(x,x') 



= [1 - /(*)] 



[ib 



i^/(x + <5 i ) 



(B3) 



and 



w[/;-r] = /(x) 



[i f + (l- s) 



(B4) 



where we sum over z neighboring unit cells at vector dis- 
tances 8i from the cell centered on x. For a cubic lattice 
Si = axi, where x^ are Cartesian unit vectors. 

To construct the master equation for the probabil- 
ity distribution functional P[f, t] of the coarse-grained 
field, we exp and in l/fi. This Kramers-Moyal expan- 
sion [3l|, l45l Hfjj leads to a Fokker-Planck equation (or 
"Kolmogorov forward equation" ) for the voter model dy- 
namics. The rates w(/;±r) will vary rapidly with the 
displacement r. However, if we choose a f2 large enough, 
we can expand oj(f] ±r) around r — in powers of 
[4^| . To second order in the master equation for 

P = P[/(x), t] follows from 

d t P = J d d x' [/ + r(x'); -r(x')] P [f + r(x'), t] 
+ W [/-r(x / );r(x')]P[/-r(x'),t] 
- [cj (/; r(x')) + u> (/; -r(x'))] P[/, t] 1 , (B5) 
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Eq. (|B6[) is a Fokker-Planck equation for the coarse- 
grained dynamics of the voter model in d-dimensions. 
For a number of applications, it is convenient to shift to 
an equivalent Langevin description Upon rescaling 
space to eliminate ft and assuming s,/xj,^h -C 1, we use 
the Ito prescription [3l[ to move from Eq. (|B6|) to the 
stochastic differential equation for /(x, t) given by Eq. © 
and Eq. f| 10[) in the main text. The Langevin equation, 
when interpreted by the rules of the Ito calculus [31 



belongs to the universality class of directed percolation 
when hi, = [20[, i.e. for "one-way" mutations to the less 
fit phase in the presence of number fluctuations. Within 
an e-expansion about d — 4 spatial dimensions, the terms 
we have ne glec ted are irrelevant in the renormalization 
group sense [20] . Note also the equivalence to a stochastic 
stepping stone model with one individual per deme [l(| ■ 
Finally, we consider the connection between our radial 
expansion models and experiments. Both the diffusion 
coefficient D r = a 2 /z r r g and the front velocity v could 
vary for different microbial radial expansions Al- 
though the Langevin description in Eq. (fT4")l and Eq. (fT5j) 
can accommodate this variation, our Bennett lattice 
model is constrained to have D r /(va r ) rj a r j (z r VT g ) < 1. 
A likely system that satisfies the bound is a Pseudomonas 
aeruginosa radial expansion on a Petri dish [l|. However, 
the bound will not be satisfied by all microbial radial ex- 
pansions. 



Appendix C: Single Sector Dynamics 

This appendix supplements the discussion in Sec. lIII 01 
by further developing the theory of the dynamics of a 
single green cell sector evolving in an otherwise red cell 
population under the radial Domany-Kinzel dynamics 
(see Sec. [H|) of a range expansion with a uniform cir- 
cular front. We will ignore mutations ([if = fib = 0) and 
allow for the green cells to enjoy a selective advantage 
s. As discussed in the main text, the Langevin equa- 
tion (Eq. ((9])) for the coarse-grained cell density /(x, t) 
is difficult to analyze for nonzero s. Hence, we move to a 
"dual" description of a single sector and consider the dy- 
namics of the two sector boundaries at angular positions 
4>i(t) and (f>2(t) at time t. Without loss of generality, let 
4>x(t) > 4>2(t) so that 4>\(t) — 4>i(t) is the angular sector 
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Suppose that the dynamics of each boundary is inde- 
pendent, and that during each generation we either gain 
or lose a single cell at each boundary due to local com- 
petition with adjacent red cells. Since the cells will gen- 
erally be staggered along the population frontier, there 
will be slight variation in the change in the sector size 
when the sector loses or gains a cell. Hence, we define a 
parameter a r equal to the average change in sector size 
when we add a single cell at the sector boundary. The 
parameter a r should be close to the effective lattice spac- 
ing a r . The master equation for the probability P{fa,t) 
of observing a boundary at fa at time t reads 



dtP(<h,t) = —P[4>i 



R(ty 

1 ~PG 



P 



R{t) 



(CI) 



where r g = 1 is the generation time, and we approximate 
[P{<f>,t + T g ) - P((f>,t)]/T g by d t P((f>,t). The probability 
of a green cell out-competing a red cell is given by pa- 

Recall that our Bennett model lattice is constructed 
so that the circular population radius advances by ap- 
proximately one cell diameter per generation. Thus, the 
angular position of the sector boundary at the frontier 
is determined by the color of approximately one cell. If 
a single green cell competes with adjacent red cells, the 
Domany-Kinzel rules (Eq. (fTJ)) imply that for s « 1, 



1 



(C2) 



where 7=1 for simple competition between a pair of 
green and red cells. The disorder in the lattice, however, 
may change 7 slightly as there can be a variable number 
of total parents ut that determine pa (see discussion 
of Eq. (fTJ) in the main text). The continuum limit of 
the master equation Eq. (|C1[) yields the Fokker-Planck 
equation, namely 



dtP(<t>i,t) = - 



jsa r 



d^P{fa,t) 



+ 



2r g [R{t)}' 



■d^P^t). (C3) 



Note that the other boundary position fait) also obeys 
Eq. (|C3|) . The Fokker-Planck equations can be converted 
to stochastic differential equations for fa (t) and fa (t) us- 
ing standard techniques [31(. Subtracting the stochastic 
differential equation for fa (t) from the one for fa (t) yields 



(C4) 



where fat) = fait) — fait) is the angular sector size 
and r)(t) is a Gaussian white noise (with (ri(t)r){t')) = 
Sit — t')). Upon defining the dimcnsionless conformal 
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FIG. 21: The speed w{s) of a green sector boundary at the 
frontier with various selection advantages s, invading an oth- 
erwise red population at the frontier. The best fit of the speed 
to a linear function for s < 0.1 {w{s) ~ 1.1s) is shown in red. 
The initial angular sector size is fa ~ tt/2, and the initial 
homeland radius Ro = 300. w(s) is calculated by averaging 
{{fat)) — fa)/ ln[R(t)/Ro] over approximately 10 3 generations, 
where {fat)) is the angular size of the green sector at time t, 
averaged over about 3 x 10 4 simulation runs. 



time r = t c {t)/t* = vt/{R + vt) = vt/R{t), we find that 
the Fokker-Planck equation associated with Eq. (|C4|) (af- 
ter changing variables from t to r) must be the one given 
by Eq. (|31[) in the main text, with 



ja r s 



and 



R vt„ 



(C5) 



This simple model is not quite faithful to the simu- 
lations on an amorphous Bennett model lattice because 
the disorder in the lattice will introduce variability in the 
constant a r and the exact number of cells at the bound- 
ary of the sector per generation. Some of these effects are 
modelled by 7, but we expect a slight s dependence in A, 
as well. Also, the cell density of the Bennett model lat- 
tice decreases a little with increasing population radius. 
However, by fitting the parameters a r and 7 to simulation 
results, we partially compensate for these complications, 
as confirmed by checks of the analytic results against our 
computer simulations. 

When s = 0, Eq. (|3ip in the main text reduces to a 



simple diffusion equation with diffusivity A 



I /Rot 



which can be solved exactly. This solution is fitted to sim- 
ulation results to find a r . We find a r m 0.91 for Ro = 300 
(recall that v = 1 in the simulations). Next, we pick a 
large initial sector angle fa ~ 7r/2 and measure the aver- 
age angular size {fat)) of the green sector as a function 
of time t. Since the sector never goes extinct in this 
case, the dynamics of {fat)) is controlled by the selec- 
tion bias. The green sector then forms, on average, the 
logarithmic spiral shape described by Eq. (|3"5j) (see also 
Ref . [r| ) . The shape of the spiral in Eq. (f3"51) implies that 
w = i(fa[t)) ~ fa)/ (In R{t) — In Ro) is a time independent 
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constant. The slope of the plot of w versus s yields an 
estimate for "fa r /v. As expected from our simple model, 



we find in Fig. [21] a linear relationship between w and s 
for s < 1 that yields 7 ps 1.2 
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